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W.H. Weatherill, F.E. Ehlers, and J.D. Sebastian 
The Boeing Company 


1.0 SUMMARY 


The problem of generating unsteady transonic air forces for use in flutter analyses remains a 
significant problem. One of the most promising procedures for solving this problem is re- 
ported by Ehlers in reference 1 . It consists of a finite difference solution for the differential 
equation of the unsteady velocity potential. The differential equation is linear and includes 
coefficients that vary with respect to space, being a function of the steady flow velocity 
potential. 

The differential equation, together with the boundary conditions, is derived in detail by 
Ehlers in reference 1 . Ehlers also presents a finite difference solution scheme based on that • 
used by Murman and Cole (ref. 2) and Krupp and Murman (ref. 3) for the solution of steady 
transonic flow. The work of reference 1 resulted in the development of the method, a pilot 
program for two-dimensional flow, and the calculation of several examples including both 
the flat-plate and NACA 64A006 airfoil. The work of this report is a direct extension of the 
earlier work and includes an investigation of solution parameters in order to reduce the com- 
puter resources needed to produce converged results; an extension of the two-dimensional 
examples of reference 1 ; the development of a pilot three-dimensional program; and an 
analysis of the dependence on frequency of the convergence of the solution scheme. 

The main results of the study are as follows: 

a) The number of iterations to solution convergence is sensitive to the value of overre- 
laxation and underrelaxation factors. 

b) There is an upper limit on frequency, depending on Mach number and size of mesh 
region, above which the relaxation procedure will not converge. 

c) Row line relaxation is more efficient than column relaxation except at combinations 
of Mach number and frequency for which convergence is marginal. 

d) The direct solution is fast and efficient for problems with a small number of grid points. 
However, the storage requirements are large and incore versions are impractical for 
realistically sized problems, even for cases in which the flow is all subsonic. The direct 
solution may provide a means for avoiding or getting around the frequency limitation 
problem discussed in b) above. 



e) A number of two-dimensional examples were calculated at Mach numbers of 0.8S and 
' 0.9 and included flat-plate and NACA 64A006 airfoils. For the flat -plate cases, the 

finite difference results compare favorably with results from linear theory using the 
program of Rowe et al. (refs. 4 and 5). Coirelation of the airfoil results with the ex- 
perimental data of Tijdeman and Schippers (ref. 6) is about the same as Ehlers foimd in 
reference 1 . In addition, two Freon calculations were made at M = 0.9. 

f) No advantage was found in using convergence acceleration methods based on the 
Aitken-Shanks delta-square process. 

g) A pilot three-dimensional program was developed for rectangular wings. Pressure dif- 
ference coefficients are presented for an aspect ratio 5 planform in harmonic pitch. 
Results for a flat-plate configuration compare well with corresponding results from 
linear theory. Results for a NACA 64A006 configuration appear reasonable, although 
there are no experimental data available for correlation purposes. 



2.0 INTRODUCTION 


An initial attempt to apply finite difference procedures in the solution of the unsteady tran- 
sonic flow problem is described by Ehlers in reference 1 . This report is a sequel to that work, 
and includes the investigation of solution parameters in order to reduce computer resources 
required to obtain useful results, an extension of the two-dimensional examples of reference 
1 , and development of a pilot three-dimensional program. 

Shortly after the publication of reference 1 , Traci et al. (ref. 7) published a paper describing 
a second solution to the unsteady transonic problem using finite differences. They, however, 
obtained a different differential equation and boundary conditions by retaining only the 
first-order time derivative term of the differential equation. The finite difference solution 
procedure appears essentially the same as that used in reference 1 . 

The purpose of this continuing investigation is to provide a practical analytical procedure 
for predicting the aerodynamic forces for flutter analyses. The examples of reference 1 
showed relatively good correlation with linear analytical results and experimental results. 
However, the amount of computer resources required to generate the results were large 
enough to make the procedure impractical for flutter calculations. This particular problem is 
significantly reduced with the use of high values of the overrelaxation factor (ORF) and 
sequential refinement of the finite difference mesh, as will be discussed in section 5.2. 

In the present work, the unsteady transonic flow is analyzed by solving for a scaled per- 
turbation velocity potential, ip. The velocity components of the flow, corresponding to the 
physical coordinates x^^ , y^ , are 

u=u^,(l+0x^) v = u„0y^ 

where u^ is the freestream velocity and 0 is the perturbation velocity potential. The scaled 
potential, ip, is related to the full potential, <j>, by the relation 


(f) = ep 


where e is assumed to be a small quantity defined in terms of the airfoil thickness ratio. 

The differential equation for the velocity potential in unsteady transonic flow as derived by 
Ehlers is 


Where 


{[K - (7+l)^o J + ^lyy 

+ ^£u^/e - "^1 = 0 


( 1 ) 


= steady velocity potential 
= unsteady velocity potential 
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M = Mach number 

CO = reduced frequency 

i =VT 

7 = specific heat ratio 

6 = (5/M)^^ where 5 is the thickness ratio or measure of camber 

and angle of attack 

K =(1-M^)/M2e 

x,y,z are the scaled coordinates x = Xo , y = HYo , z = Mz© > where Xo , y© > Zq 
are the physical coordinates and p = 5*''^ M^''^ 


Equation (1) may be rewritten as 

fu'P, - (2ioo/e)</> 1 +<P. +'P. +q'^j = 0 (2) 

L * X U X * yy * zz * ^ ^ 

and this will be the form used throughout the report. For two-dimensional flows, equation 
(2) is rewritten as 


- (2ico/e)^j]^ +'^iyy+q'^i = 0 

The wing shape as a function of time is written 

Zq = 5f(x,y,t) 


(3) 

(4) 


and the linearized boundary condition for the total velocity potential, 'P , is 

<^z = fx<^x,y,t) + ft(x,y,t) (5) 

For two-dimensional flow, the is replaced by y^ and f = f(x,t). 

Equations (2) and (3) are rewritten as finite difference equations and solutions are obtained 
using relaxation procedures. The finite difference equations together with far-field boundary 
conditions for three-dimensional flow are presented in appendices A and B. 

An important development during the course of the program was the encountering of dif- 
ficulties in obtaining solution convergence-first for row relaxation, and then for certain 
combinations of Mach number and frequency. This led to analytical investigations of the 
relaxation procedures used that are discussed in section 4.0 and described in detail in ap- 
pendices C through E. 
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Empirical investigations of the relaxation procedures are presented in section 5.0 with two 
dimensional examples presented in section 6.0. 

Finally, section 7.0 describes the three-dimensional pilot program together with some 
examples. 


5 




3.0 SYMBOLS AND ABBREVIATIONS 


a,b 

ADI 

b 

BSOR 

c,d 

c 

Cs, , Cgj , dj, , djj 

Cl , di , Cj , dj 
ACp 

E 

ERROR 

f(x,y,t) 

f(x,t) 

fo 

f, 

Fij 

h 

ij\k 


coefficients for y, z differences corresponding to seooiyl derivatives, 
with appropriate subscripts (eq. A3); also length of sides of region 
for finite difference solution 

altemating-direction-implicit iteration scheme 

semichord of wing 

block successive overrelaxation 

coefficients for x difference corresponding to second derivative 
(eq. A3) 

height of side of region for finite difference 
equation (A1 3) in appendix A 
equation (A27) in appendix A 

coefficients for second-order accurate difference corresponding to 
first derivative (eq. A3) 

jump in pressure coefficient 

coefficients in difference equations with appropriate subsaipts; also 
used as unknown error in Von Neumann analysis 

see equation (19) 

instantaneous wing shape defined by Zo = 5 f(x,y,t) 
instantaneous airfoil shape defined by yo =5 f(x,t) 
undisturbed wing or airfoil shape 
unsteady contribution to wing or airfoil shape 
see equation (A6) 

^k_ + 1 * ^k„ 
m m 

x,y,z subscripts for points in the mesh 




X index for first mesh point behind hinge 

imax>iihax» ^^max 

maximum number of x,y,z mesh planes respectively 

ii 

X index for trailing edge 

Iw 

integral defined by equation (B2) 

im 

y mesh line just below airfoil (two-dimensional configuration) 

K 

(1 -M2)/M^e 


z mesh plane below wing (three-dimensional configuration) 

LfJ 

superscripts denoting lower or upper boundary on Fy ; also used to 
denote lower and upper triangular matrices resulting from matrix de- 
composition (sec. 5.3.2) 

m 

subscript used on mesh point indices to denote points adjacent to 
and below the airfoil or wing 

M 

freestream Mach number 

n 

number of iteration 

Nx,Ny,Nz 

number of mesh increments in x,y,z directions 

ORF 

overrelaxation factor 

P(x,y) 

ip, + ioxpj acceleration (or pressure) potential 

q 

- i<o(y- O'Poxx 

r . 

overrelaxation factor 

R 

variable defined for equations (A1 8) 

Ri , Ro 

variables defined in equation (B4) 

s 

semispan of wing 

Si 

<>Sn, + 2 - >5m 

Sj 

<>'im - S'im-l*''' 

u 

R - (-rH)«>ox 

URF 

underrelaxation factor 
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Xo,yo>Zo 


physical coordinates 


x,y,z 


x'l, y\, z'l 


x,y,z 

X£, x^ 

*a 

Ax, 

Ax, 

Vt 

/3 

Pi . P2 
7 

A0I 

e 

S 

62 

X, 

/i 


X 


scaled coordinates (xo, Myo> M^o) for the three*dimensk>nal problem; 
the scaled coordinates for the two-dimensional problem ai« x and y, 
with X again being the direction of fluid flow 

' J 

variables of integration 
coordinates defined in equation (A1 8) 
coordinates of leading and trailing edges 
coordinate of control surface hingeline 

Xj-Xj., 

*i-i - *i-2 
coordinate of wing tip 

n/1 - 

parameters in equation (Cl 8) 

ratio of specific heats for air 

jump in 0] at plane of wing or vortex wake 

( 5 /M) 2/3 

thickness ratio or measure of camber and angle of attack 

X, - X, 

*i ~ *i -1 
‘max 'max ' 

toM/(l - M^) 

Scale factor on y® and Zq, = . i 

unsealed perturbation velocity potential 

steady scaled perturbation velocity potential. 

unsteady scaled perturbation velocity potential 

wake integral defined in equation (A24) 

acceleration or pressure potential 
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fundamental source solution of integral equation for evaluation of 
far-field boundary conditions 

angular reduced frequency (semichord times frequency divided by the 
the freestream velocity) 

critical reduced frequency (sec. 4.2) 

Kronecker product of matrices A and B (see sec. D.l .3 in app. D) 



I 


4.0 ANALYTICAL EVALUATION OF THE CONVERGENCE 
OF THE RELAXATION TECHNIQUES USED FOR 
THE UNSTEADY TRANSONIC FLOW PROBLEM 

During the course of the program, a considerable effort has been made to analyze and im- 
prove upon the relaxation procedures used to solve the finite difference equations. At 
first this was done to ensure the convergence of the mixed-flow problem using row relax- 
ation. It was found that a straightforward implementation of this scheme, in direct analogy 
to the column relaxation scheme, was divergent for the mixed-flow problem. Following the 
recommendations of Jameson (refs. 8 and 9), an analysis was made in which successive itera- 
tions were treated as a pseudotime variation. This analysis showed the need for the addition 
of time-like difference terms in the program. A summary of this work is given in section 4.1 
and a detailed discussion is presented in appendix C. 

The next problem was encountered in attempting to complete some of the numerical 
examples. In particular it was found that, for a given Mach number, there was an upper limit 
on the reduced frequency at which the relaxation solution converged. At higher values of 
frequency, the solution was found to diverge. This behavior was found to be essentially 
independent of the type of relaxation scheme (i.e., row or column), the mesh spacings, the 
relaxation factor, and whether the differential equation was purely elliptic (flat-plate con- 
figuration) or of mixed elliptic and hyperbolic type (airfoil section configuration). A matrix 
analysis of a simplified version of the problem revealed that there is indeed an expected 
frequency limitation on convergence, and this may be recognized by noting the similarity 
between the unsteady differential equation and the Helmholtz (or reduced wave) equation. 
This analysis is summarized in section 4.2 and a detailed account presented in appendix D. . 

Finally, a third means of analyzing relaxation procedures is the Von Neumann stability 
criterion. The application of this analysis method to the pertinent relaxation process gen- 
erally concurs with the findings of the preceding two analyses. The details of this analysis 
are presented in appendix E and summarized in section 4.3. 

4.1 A TIME-LIKE CHARACTERISTICS ANALYSIS 


Jameson treats the difference between two consecutively iterated values of the velocity 
potential as a time derivative, namely: 


^(n -1) 


- = -At 

ij 


g,p(n) 

*ij 

at 


( 6 ) 


where the superscripts denote the iteration number and the subscript i, j denotes the finite 
difference point. When this expression is substituted into the difference equation and all 
terms are expanded in a Taylor’s series about the central point, a differential equation 
results that contains not only the terms of the original differential equation but also addi- 
tional time derivatives. A study of variation with time of the solutions to this differential 
equation with arbitrary initial conditions will reveal some insight into how the iteration 
method will converge. 
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For flows that are subcritical, the simple straightforward difference equation developed by 
Murman and Cole (ref. 2) will converge whether column solution or row solution of the 
difference equations are used in the relaxation procedure. Column relaxation will also work 
for mixed flow, but row relaxation will diverge unless time-like difference terms are added. 
The time-like differential equations for both column and row relaxation are derived in 
appendix C. At elliptic points, row relaxation for sweeping toward the airfoil (y < 0) leads 
to the differential equation (C 1 5), i.e. 


(u^ 


) .JM(P +ip +qiPj +At • [23: --^(a.+ b.)]'^. 

‘ X X e * X * yy * L J r J J J ‘ t 


(7) 


where Ayj = yj + 1 - yj - ] and r is the overrelaxation parameter. 'Fhe differential 

equation is expressed in canonical form by introducing a new time t = t + — y. 

Ayj 

Thus 

')\i \ r2(a;j.b:) -| 

+ 'Pi + q<P, - At • 2a: <P 

*xx e*x *yy ‘ Lr JJ‘r 


(w) 




= 0 


"TT 


( 8 ) 


This equation is hyperbolic in the time and hence initial value problems are properly posed. 
For supersonic points, u < 0 and terms have the same sign rendering the equation 
hyperbolic in y. Since supersonic points experience only upstream influences, the equation 
must be made hyperbolic with respect to the variable x. This is achieved by adding differ- 
ence terms of the form yielding time derivatives and ^ and choosing their coef- 
ficients so that the resulting coefficient of the v^j term is positive. The method was em- 
ployed by Jameson for the steady-state equation for the velocity potential and was found to 
be successful for the linear unsteady differential equation for harmonic motion discussed in 
this document. The inclusion of these additional terms does not affect the final solution 
since these terms become negligible when the solution is converged. 


Introducing the time derivative terms into the difference equations for column relaxation 
at elliptic points leads to the following time-like differential equation for sweeping in the 
increasing x direction: 

{u<p, ) -•^^'Pi + 'Pi - 2 ua 'P - 2 (ua -iojoc )'P = 0 

i^x € *x \yy ^ ‘ ^xt ^ ^ 't 


where 0:1,0:25 arid 0:3 are given in equation (C8). With the new time variable, r = t + o:^ x, 
the equation takes the canonical form 


(UV’j ) 
*X X 


2ioo 


P, + P, 

* X 




ua, pi, =0 


( 10 ) 


TT 
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Since for elliptic or subsonic points, u > 0, this equation is hyperbolic in time. At hyper- 
bolic points, backward differences are used for the x derivatives and the resulting differen- 
tial equation contains no 'Pi term but takes the form (eq. Cl 1 ) 


(uV?. ) - +<P +q<^, -2u-jr — 

* X X e * X ' yy * AXj (AXj + AXj )\ t / > t 

+2i^fLAt i-. V. =0 

e \ Ax, . Ax, + AXj / 4 


where Axi = xj - Xj . i and Ax 2 = xj . j - xj . 2- 


( 11 ) 


Since u is negative, this equation is parabolic in time, but hyperbolic with respect to x. To 
have damping at supersonic points, equation (1 1) indicates that r must be less than 1 , and 
hence underrelaxation must be used at supersonic points. 


4.2 A MATRIX CONVERGENCE ANALYSIS 


The observed frequency-dependent limitation on the convergence of the overrelaxation 
method may be analyzed according to a system matrix approach. In outline, the method 
proceeds as follows: first, simplifying assumptions are made that replace the original 
problem with that of solving 




XX 


K 


yy 


2ico 

eK 




CO 


eK 


= 0 


( 12 ) 


over a rectangle of sides a, c and with being prescribed on the sides. Then the region and 
differential equation are discretized using a uniform mesh and central differences, respect- 
ively, leading tq_^a systernpf linear difference equations. This system is written in matrix 
form as A^i = R, where is the vector of unknown values of at the interior mesh 
points and R is a vector containing the boundary values. The elements of A are known 
functions of the reduced frequency, co. A theorem is invoked that states that under certain 
mild assumptions, the line or block overrelaxation scheme used in the program will converge 
if, and only if, all the eigenvalues of A are positive; i.e., A is positive definite. To apply the 
theorem, the eigenvalues of A, which are, of course, also functions of oo, are determined; 
then the value of co for which the smallest eigenvalue becomes zero is found. We call this 
value the “critical frequency” denoted by cocr, since by virtue of the theorem, it is the 
value of CO below which relaxation will converge and above which it will diverge. 


Exact and approximate formulas for co^r derived in this way are given in appendix D. Here 

we give only the approximate formula, which is 

2 - _ 

1-M 

M 


CO 


cr 


= 7T 


[a Kc J 


(13) 


where a and c are the width and height of the mesh region. This completes the outline of the 
analysis, the goal of which was to find a formula for coct to explain the frequency limita- 
tion. We now turn to validation of the formula. 
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A comparison of the critical frequency predicted from the analysis of the simplified prob- 
lem with various computational observations is given in figure 1 . The solution region in all 
cases was defined as the rectangle x = -2.65 to 2.75; y = -6.25 to 6.25. The approximate 
prediction formula given by equation (13) was used since the number of mesh points em- 
ployed (at least 25 x 20) was sufficiently large in each direction so that the difference be- 
tween ojct from the exact and approximate formulas is less than 1%. 

With regard to the computational observation points, those indicating convergence corres- 
pond to the largest frequency for which convergence was obtained at the given Mach num- 
ber. In the two airfoil cases, the convergence was marginal and resulted only after substan- 
tial experimentation with the values of the under- and overrelaxation factors, at supersonic 
and subsonic points, respectively. In the case of the flat plate, additional computations were 
performed at slightly higher frequencies; actual divergence of error measure was observed, 
as indicated in figure 1 . 

Further validation of the results of the matrix analysis was obtained as follows. It is an im- 
mediate inference from equation (13) that a decrease in the dimensions of the solution 
region implies an increase in the frequency for which convergence of overrelaxation can be 
obtained. As a test of this inference, and hence of equation ( 1 3), the dimensions of the 
solution region were reduced from 5.4 to 2.8 in the x -direction and 1 2.5 to 6.0 in the z- 
direction. For these dimensions and M = 0.9, the critical frequency predicted by equation 
(13) is 0.254. Flat-plate computations were perfonned for co = 0.25 and co = 0.30, with 
the result that convergence was observed in the former case and divergence occurred in the 
latter case. 

The existence and general location of a critical frequency, dependent on Mach number and 
predicted by the analysis of the simplified problem, are in good agreement with the com- 
putational results from the full problem. Prediction of the exact location of the critical 
frequency is not to be expected since, of the five assumptions made in the formulation of 
the simplified problem in appendix D, all in the airfoil case and b through d in the flatplate 
case are violated in the actual computational scheme. The results strongly indicate that 
the cause of the frequency limitation in the full problem is the same as that in the simplified 
problem: the failure of the system matrix to remain positive definite. 

It should be noted that the conclusions of this section so far are relevant to relaxation 
solutions. For a direct solution where A is formed from the complete set of simultaneous 
albegraic equations and the solution is obtained by matrix inversion, the matrix A need only 
be nonsingular rather than positive definite for solutions to exist. Thus use of a direct solu- 
tion procedure may well lead to solutions at reduced frequencies above oocr- However, since 
the far-field must be updated as the velocity potential distribution changes, the matrix form 
is more nearly A.(p\ = R«^i ‘ \ where and ‘ H are the vectors of values of 

the unknown velocity potential and the values determined in the preceding solution, 
respectively. This iterative form resulting from the calculation of the far-field boundary 
conditions imposes the additional condition for solution convergence that the effective 
eigenvalues of the matrix product A"' R must be less than one in modulus. 
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4.3 A VON NEUMANN STABILITY ANALYSIS 


A third means of examining the convergence of a differencing method is the Von Neumann 
analysis of error propagation. Let be the error at the k, £ grid point; then for a uniform 
X and a uniform y mesh the errors may be expressed in the form 

P ‘*^max ^ max 

where a and c are the width and height of the mesh region. Since the difference equation 
is linear, only a single term need be analyzed. 


A solution of the difference equation of the form 


{p (n) 

‘k£ 




where n is the order of iteration, has the appropriate initial value consistent with equation 
(14). In order for errors to decrease with increasing n, g must have a magnitude less than 
one. In appendix E, the equations of column relaxation, row relaxation, and an ADI (alter- 
nating direction implicit) method were analyzed. The condition I gl < I in all cases reduced to 
the same inequality limiting the range of frequency, namely: 
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where d, = 27rp^^and0j = 2irq'^. 

riiis relation sliows that the range of frequency is increasingly limited as the Mach number 
approaches 1 . 


In the hopes of obtaining better convergence for the higher reduced frequency problem, an 
ADI (alternating direction implicit) method was tried. In the flat-plate case, we consider for 
the timcHlependent equation for two-dimensions 
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A uniform mesh is swept through using a row solution of the difference equation 
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and then through the mesh by a column solution of the difference equation 
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Applying the Von Neumann stability analysis to the preceding two equations yields the 
same inequalities as found for conventional row and column relaxation. This indicates that 
the ADI method yields no improvement in the convergence for higher frequencies. The ADI 
method was coded and tried but failed to converge. Lack of time has prevented an adequate 
investigation of this failure. 
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5.0 SOLUTION PROCEDURE INVESTIGATION 


An empirical investigation was conducted of the finite difference solution procedure to 
determine the effect of va.rying solution parameters and techniques and in minimizing solu- 
tion computer resources. The investigation included: 

• Relaxation factors 

• Grid point distribution and spacing 

• Extent of mesh area 

• Sequential refinement with respect to Mach number, reduced frequency, and number 
of mesh points 

• Row relaxation 

• Direct solution 

• Convergence acceleration 

The results of these investigations are summarized in this section and detailed in appendix F. 

5.1 SOLUTION CONVERGENCE CRITERIA 

For this report, solution convergence was determined by monitoring the ERROR, which is 
defined as maximum value of all i, j, k of 


(n) _ ^(n-1) 
* jjk ‘ jjk 


(19) 


where i^| is the unsteady velocity potential for the nth iteration, is the corres- 

|X)nding potential for the preceding iteration, and r is the relaxation factor. The solution 
was considered converged when ERROR = 10‘® . In some cases, particularly for finer 
meshes and pitch mode, convergence was considered complete when ERROR < 10'^. 


The maximum residual was considered to have greater potential as an indication of solution 
convergence than ERROR. The residual is a measure of the degree to which a solution (here, 
a set of velocity potentials) satisfies the finite difference equations. The iteration may be 

written in general matrix form as A^j = R^j . The residual at the ith point is thus 
the summation £(Ajj - Rij)'^ij- Preliminary attempts to calculate the residual resulted in 


values several orders of magnitude larger than the corresponding ERROR. Multiplying 
through by an area associated with each mesh point reduced the difference between the two 
values to some two orders of magnitude. However, the precise relationship between ERROR 
and the maximum residual is not yet entirely clear and thus further investigation is required 
before the residual can be used as a convergence criterion. 
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In the discussion on relaxation factors it will be noted the value of ERROR may be reduced 
rapidly by reducing the value of r. However, the resulting pressures were inconsistent with 
the minimal level of ERROR that had been reached with a single value of r. This again indicated 
that some form of residual measure might provide a more realistic insight into convergence. 

5.2 SOLUTION PARAMETERS 
5.2.1 RELAXATION FACTORS 


In the iteration solution procedure used for this report a relaxation factor was applied to 
the velocity potential using the following replacement formula 





(1 -r) 


Mjk 


( 20 ) 


where is the velocity potential form the previous iteration at point (ij,k), and 

on the right-hand side is the result of the current iteration. The factor r is the relaxation 


factor and is set to some value between 0 and 2. The procedure here was to use overrelaxation 
( 1 .0 < r < 2.0) for points at which the steady flow was subsonic and the unsteady differential 
equation was elliptic, and to use underrelaxation (0 < r < 1 ) for points at which the steady 
flow was supersonic and the differential equation was hyperbolic. 


The characteristics of convergence for overrelaxation will be discussed with figures 2 and 3. 
These sketches show ERROR versus number of iterations. It is convenient to plot the ERROR 
in terms of log scale and the number of iterations as a conventional linear scale. It is noted that 
a straight difference of successive potential values includes the relaxation factor as a multiplier. 
Thus, we divided r out in calculating ERROR to provide a “normalized” measure of the error. 
I'or small ORF’s (overrelaxation factors) these curves were made up of two nearly linear 
portions; for the initial iterations the convergence was rapid, resulting in a steep slope during 
the early iterations and a very shallow slope for the later iterations. As the ORE was raised, 
the .slope.of the initial iterations became less, while that of the later iterations increased. This 
trend continued as ORF was raised until the convergence during the initial iterations eventually 
became mildly unstable, but still converged well. Finally, when ORF was too high, convergence 
became ()uite unstable, and the slope of ERROR versus number of iterations decreased sharply. 

An ORF should be selected, if possible, so that the convergence criteria is attained within 
the number of iterations for which the steep initial slope holds. 


Another interesting characteristic (fig. 3) is that if the ORF is changed during a solution 
calculation, there is an immediate change in the ERROR versus number of iterations curve. 
If the ORF is lowered, there is a significant drop in the ERROR curve. However, after the 
drop, the line levels out at a slope less than the slope for the higher ORF. It is possible, of 
course, to take advantage of this phenomenon to achieve a specific ERROR value for which 
convergence has been defined. However, we encountered inconsistent results where the con- 
vergence ERROR had been set as low as 10'^ , and thus caution is advised in making use 
of this characteristic. 
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The selection of the optimum overrelaxation factor was dependent on: 

(a) The number of mesh points 

(b) Level at which the convergence criteria was set 

(c) Relaxation procedure (row, column, etc.) 

(d) Solution sequence (start at upstream boundary, trailing edge, etc.) 

(e) Frequency for a given Mach number 

Values of the relaxation factor ranged from 1 .9 down to 1 .4 depending on the particular 
example. Item (c), the variation with frequency, appeared to be involved with the frequency 
limitation problem and is not discussed in appendix F. However, for a given Mach number 
and as the frequency was increased into the range of marginal stability, the. value of the op- 
timum relaxation factor decreased. 

Generally, the coarser the mesh, the longer the initial slope lasted, and thus the smaller the 
optimum ORF. In this sense, the optimum ORF was dependent on the level at which the 
convergence criteria was set. 

The pattern of the variation of optimum ORF with relaxation procedure and solution 
sequence was not explicitly defined within the limited number of examples run. Some 
illustrations of the variations encountered are presented in section F.l .1 of the appendix. 

The use of complex ORF’s does not improve the solution convergence characteristics. This 
was first demonstrated empirically and later shown with the Von Neuman stability analysis 
presented in appendix E. 

The selection of ORF is best done by trial and error since it is a function of so many vari- 
ables. We would recommend as large an ORF as possible and we have had considerable suc- 
cess with ORF’s of 1 .85 and 1 .9, as shown in the appendix. 

The solution convergence was not generally as sensitive to the underrelaxation factor as to 
the overrelaxation factor, although one case is cited in appendix F when the solution 
diverged for an URF of 1 .0 and converged rapidly for an URF of 0.7. 

5.2.2 GRID DISTRIBUTION AND SPACING 

The examples of grid distribution and spacing showed that the representation of the 
pressures in the neighborhood of the flow singularities was significantly improved by cluster- 
ing the points about the singularities and increasing the number of grid points. 

For the first case, the mesh-area dimensions were fixed and the total number of points held 
nearly constant. The mesh points were spaced such that the ratio of sizes of adjacent intervals 
was X, where X could be varied between 2/3 and 3/2. The smallest intervals (for X > 1.0) 
were centered about the known flow singularities; i.e. at the wing leading edge and the con- 
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trol surfiRL* hingeline in the flow direction and at the airfoil centerline in the crossflow 
direction. Successively improved representations of the leading-edge singularity in the 
pressure distribution was obtained as A changed from 1 .05 to 1 .25 to 1 .4. At the hingeline, 
tlic pressure distribution changed significantly between X = 1 .05 and 1 .25 and minimally 
between X = 1 .25 and 1 .4. 

For the second case, the mesh area dimensions and the point spacing factor, X, were held 
fixed, and the number of mesh points varied. A significant improvement in the pressure 
representations was obtained in going from a mesh of 25 x 1 6 to a mesh of 34 x 28. The 
change in pressure distribution in going from a 34 x 28 to a 42 x 30 mesh was minimal. 

5 .2.3 EXTENT OF MESH 

The effect of varying the extent of the finite difference mesh was investigated by altering 
the location of the upper and lower mesh boundaries. Solutions were obtained for the 
boundaries at yoniax “i9,±18.5, and ±29.6 (in physical coordinates) while the number 
of mesh points was held fixed. The pressure distributions and the velocity potential were 
compared. There was surprisingly little variation in the pressure distributions for the three 
cases. However, the distribution foryoj^^x" slightly smaller in amplitude than the 

other two, which in turn were essentially the same. The velocity potential showed much 
more difference between the three cases, a difference not reflected in the pressure distribu- 
tions. The best representation of the leading edge singularity in pressure was with Yo ~ 
+9, the calculation apparently benefiting from the compression of points in the crossflow 
direction. 

5.2.4 SEQUENTIAL REFINEMENT 

In starting a new analysis, a set of zeros is often used as the initial values for the velocity 
potential distribution. The question is whether solution convergence for the desired set of 
parameters may be most economically achieved by calculating the velocity potential for 
intermediate values of the parameters, and using the resulting potentials as the initial dis- 
tributions for the final calculation. This process of sequential refinement was applied in 
terms of Mach number, reduced frequency, and the number of mesh points. 

The example presented in section F.1.4 of the appendix shows relatively little difference in 
the number of iterations to convergence, whether the initial velocity potential distribution 
is set to zero or taken from a corresponding problem with a relative small difference in Mach 
number or reduced frequency. The conclusion is that if potential distributions exist for 
intermediate values of Mach number and frequency, it would be worthwhile to use them. 
However, it would not be worthwhile to calculate them as an intermediate step in calculat- 
ing the potential distributions for the desired values of Mach number and reduced 
frequency. 

Sequential refinement with respect to the number of mesh points does appear worthwhile. 
Here, the velocity potential distribution is calculated for a relatively coarse grid. Then, the 
potential distribution is interpolated to a finer mesh and the resulting distribution used as 
initial values for another iteration sequence. An example is given in section F.1.4 of the 
ap|)ciHlix where the number of iterations for the refined grid is cut in half using the results 
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from a coarse grid. The actual savings is a function of the relative number of points of the 
coarse mesh with respect to the fine mesh and the number of iterations required to obtain 
a converged solution for the coarse mesh relative to the number required for the fine mesh, 

5.3 SOLUTION PROCESS 


5.3.1 ROW RELAXATION 

The finite difference equations as presented in reference 1 were written for column relaxa- 
tion where the solution is for a line of points parallel to the z -axis in three-dimensional 
flow (or the y -axis in two-dimensional flow) and extending from the lower boundary to the 
upper boundary. By rearranging terms, the equations may be rewritten for row relaxation 
where the solution is for a line of points extending from the upstream boundary to the 
downstream boundary. For subsonic flow, the resulting row formulation provides relatively 
rapid solution convergence. However, for mixed flow, additional terms must be included 
in the finite difference equation to obtain convergence. These terms, resulting from a time- 
like analysis of the finite difference equation, are derived in appendix C. The application of 
these terms is discussed in section F.2.1 of appendix F. 

Examples of solution convergence using both row and column procedures, and using several 
solution sequences (i.e., the order in which the rows and columns are solved) are also presented 
in section F.2.1. Generally, row relaxation was found to be significantly more efficient than 
column relaxation for both subsonic and mixed flow. The only exceptions to this were for 
values of Mach number and frequency for which relaxation solutions were marginally stable. 
Here, although neither procedure was particularly rapid, column relaxation provided solution 
convergence for values of reduced frequency at which the row procedure had started to 
diverge. 

Row relaxation was implemented in three ways. The most efficient manner was to start at 
the upper and lower boundaries and work toward the wing surface, alternately taking a row 
from the top section and a row from the bottom section. The alternatives of starting at the 
wing and working out toward the upper and lower boundaries and of starting at the lower 
boundary and simply taking successive rows in working to the top boundary were also tried. 
Column relaxation was run in two ways. The most efficient procedure proved to be the 
sequence that started at the trailing edge and worked toward the upstream boundary, then 
moved to the column just aft the trailing edge and worked to the downstream boundary. 

The alternate sequence started at the upstream boundary and moved by successive columns 
to the downstream boundary. 

5.3.2 DIRECT SOLUTION 

A version of the pilot two-dimensional program was written to provide a direct solution for 
the interior velocity potential distribution. By direct solution, we mean that the complete 
set of equations is solved all at once rather than in subsets, as with the row or column relax- 
tion. For the direct solution as implemented here, there is stiU an iterative loop since the 
boundary conditions are calculated using an existing set of velocity potential distributions 
(i.e., the velocity potential distribution calculated in the preceding pass). The program for 
the direct solution was set up only for the purely subsonic flow. 
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The direct solution program has been used on the standard test problem of a flat plate with 
oscillating quarter-chord control surface for grids of 1 7 x 1 0 and 25 x 1 6 mesh points. It has 
also been used for the flat-plate pitch case. The results compare almost exactly with the re- 
laxation solutions. The number of iterations required to reach convergence was less than 20 
for all cases and computing time appears minimal. 

In general, the direct solution was found to be a relatively fast and efficient solution pro- 
cedure for small problems (coarse grids), but storage requirements make it impractical for 
incore solution of two-dimensional problems with what we currently feel to be realistic 
sized finite difference grids. The upper limit would be of the order of some 700 mesh 
points. As noted in a preceding section, current studies indicate that practical grids are of 
the order of 1200 or more points. Also, storage requirements for the mixed flow problem 
are even larger. Because of its speed, however, the direct solution has provided a very useful 
experimental tool for testing modifications to the program. Use of an out-of-core direct 
solution method has not been attempted, but is one possible direction for further investiga- 
tion, particularly in light of the frequency limitation problem discussed in section 4.2. 

5.3.3 CONVERGENCE ACCELERATION METHODS 

The relatively regular, uniform, and monotonic behavior of the pressure difference distribu- 
tions with successive iterations has suggested the use of convergence acceleration techniques. 
Also, these procedures have been successfully applied in limited examples of steady transonic 
flows by Hafez and Cheng (ref. 10) and Martin and Lomax (ref. 1 1 ). Further, our studies 
have shown relatively good behavior of both the velocity potential distribution and the 
ERROR with successive iterations. Despite the optimism with which this study was ap- 
proached, the results were not favorable. 

The Aitkin-Shanks nonlinear transformation (5^ -process) was tried first and applied to the 
velocity potential distribution. Typically, the extrapolated velocity potential had a large 
value of ERROR. Additional iterations resulted in a rapid drop in ERROR back to the 
convergence path that the solution was following before extrapolation. Although little 
was lost by using the extrapolation process, nothing was gained either. 

Since examination of the extrapolated velocity potential showed it to be inconsistent with 
that which would be estimated by eye, a modified form of the Aitkin-Shanks transformation 
equation was introduced in order to constrain the shape of the extrapolation. This was tried 
on both the velocity potential and the pressure distributions without showing an improve- 
ment in solution convergence over straight relaxation. 

Generally, the results of these studies have been discouraging despite the fact that con- 
vergence appears to be monotonic for the pressure and velocity potential distribution as 
well as ERROR. Results have not been improved by working with solutions with the 
smoother convergence characteristics obtained by either reducing the overrelaxation factor 
or using results from higher numbers of iterations. Finally, the real part of the solution ap- 
pears to behave much better than the imaginary part, indicating that convergence accelera- 
tion procedures may well be much more promising for the steady-flow problem. 
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6.0 TWO-DIMENSIONAL EXAMPLES 


A series of two-dimensional examples has been computed in order to further explore the 
accuracy and characteristics of the finite difference method examined in this report. The 
most important result from these examples was the discovery that solution convergence is 
limited as a function of Mach number and frequency. This particular point is discussed in 
section 4.0. This problem has reduced the number of examples in terms of number of dif- 
ferent frequencies examined. The examples included both the NACA 64A006 airfoil and a 
flat plate (i.e., a section of vanishing thickness). The range of Mach numbers was 0.85 to 
0.90 and the oscillatory motion included both section pitch and quarter-chord control 
surface rotation. 

Generally, the rate of convergence was very dependent on the value of the relaxation factor. 
This was particularly true for elliptic points where an overrelaxation factor was used. Indeed, 
the selection of the overrelaxation factor is significant enough so that it ought to be deter- 
mined for each case separately. Mixed-flow solutions appeared relatively less sensitive to the 
underrelaxation factor used with hyperbolic points. However, the case M = 0.9 and co = 0.06 
for the airfoil in pitch diverged with an URF of 1 .0, whereas 0.7 led to rapid convergence. 

The sequential refinement in terms of mesh size (the number of points rather than the total 
mesh area that was held constant) proved worthwhile, whereas refinement in terms of fre- 
quency did not. The solutions presented here were obtained using some 1 50 to 200 itera- 
tions of a 25 X 20 grid followed by some 200 iterations of a 42 x 30 grid. 

It was during the calculations for these examples that the limitations on convergence in 
terms of frequency and Mach number were encountered. Since the convergence problem 
existed for flat plate as well as the airfoil section, the handling of the shock and the attend- 
ant mixed flow was not the cause. At first it was assumed that the difficulties were due to 
poor selection of ORF or mesh size, or to some other parameter of the solution procedure. 
Thus considerable experimentation was done in an effort to obtain a significant improve- 
ment in the solution convergence. This included, besides ORF and mesh-size variations, 
various forms of column as well as row relaxation, the inclusion of a second-order approxi- 
mation to the far field, and variation in mesh-point spacing. Sequential refinement in terms 
of frequency did not help. In addition, a considerable reshuffling of terms from one side to 
the other in tifie finite difference equation (eq. A1 ) was tried. Also, the program was modi- 
fied to solve the differential equation and boundary conditions used by Traci et al. in ref- 
erence 7. None of these changes provided the significant improvement that we felt was 
necessary to the problem. 

As an example, consider the case of pitch for a flat plate with a harmonically oscillating 
quarter-chord control surface at M = 0.9. For a 25 x 20 mesh, relatively rapid convergence 
was obtained at frequencies up through co = 0.1 2. Although difficulty was encountered, 
convergence was obtained at co = 0.14 with column relaxation working considerably better 
than row relaxation. However, at co =0.1 6, solutions converged rapidly to an ERROR of 
10“® and then the solution curve flattens out to provide, at best, very slow convergence. 

At o) = 0.1 8, the solution diverged (this for a 42 x 30 mesh; ORF = 1 .0 under row relaxa- 
tion). By dropping theq<Pj , :term of equation (3) and simplifying the boundary conditions. 
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the program was readily modified to approximate the solution of Traci et al. (ref. 7). Here 
the solution converged for co = 0.1 6, and at cj = 0.1 8 the convergence became very slow 
again. Finally reduction of the mesh area, as suggested by the analyses discussed in section 
4.0, provides rapid solution convergence at co = 0.1 8. 

The flat-plate results are compared directly with solutions from the NASA linear subsonic 
three-dimensional unsteady aerodynamic program (refs. 4 and 5). This program solves the 
pressure-downwash integral equation using assumed pressure modes. Since the steady 
transonic flow past a flat plate is uniform, the finite difference results should match the 
subsonic results directly. Figures 4 through 8 present the results for a flat plate (section of 
vanishing thickness) in pitch. In all cases the correlation between the finite difference 
program and the NASA program is good.The finite difference method appears to underes- 
timate slightly the amplitude as calculated with the NASA program. Also, the degree of 
correlation appears independent of Mach number and frequency. Solution convergence was 
considered to be a maximum error of 10”'* between velocity potentials for successive 
iterations. This appears to be adequate for the calculations involving the pitch mode. 

Figures 9 and 10 show the jump in pressure coefficient across a flat plate with an oscil- 
lating quarter-chord trailing-edge control surface. Results are presented for M = 0.9 and 
reduced frequencies of 0.06 and 0.1 2. Again correlation between the finite difference pro- 
gram and the NASA program appears good. 

Calculations for the NACA 64A006 airfoil section are presented in figures 1 1 through 1 7. 
Results for pitch motion at M = 0.85 and reduced frequency of 0.06 and 0.24 are shown in 
figures 1 1 and 1 2. In figure 1 1 , results for three different mesh grids are shown while the 
total mesh area is held constant. First, it is noted that increasing the number of points 
(perpendicular to the flow) improves the representation of the leading-edge singularity. 
Secondly, the clustering of points about the shock, which in this case is just forward of 
midchord, results in a kind of singularity, with the pressure going to large positive numbers 
in front of the shock and reappearing from negative numbers behind the shock. For this 
example (in the 45-point distribution), points are clustered about the 3/4 chord, while in 
the 46-p6int distribution, the points are clustered about x = -0.1 , the approximate location 
of the shock center. This latter phenomenon is illustrated in figure 1 2, which presents the 
results for two flow-wise point distributions. There are no experimental or other analytical 
data for correlation with these pitch motion results. 

A two-dimensional example was also calculated to compare the change in pressure distribu- 
tion that results from using the specific heat ratio of Freon (y = 1 .1 35) instead of that for 
air (7 = 1 .4). Figure 1 3 shows the steady-state pressure coefficient at M = 0.9 for Freon and 
air for the NACA 64A006 airfoil. The shock in Freon is slightly ahead of the shock in air. 
Figure 14 shows the jump in pressure coefficient across the airfoil due to the harmonic 
pitch of the airfoil section with a reduced frequency of 0.06. 

The remaining figures (figs. 15 through 17) show pressure difference coefficient distribu- 
tions for oscillatory quarter-chord control surface motion. The cases for figures 1 5 and 1 6 
are calculated at M = 0.875, and M = 0.9 for co = 0.06, and are shown with the correspond- 
ing measured data from Tijdeman and Schippers (ref. 6). Solutions at M = 0.9 and co = 0.1 2 
were marginally convergent and converged pressure distributions were not obtained. Correla- 
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Figure 7. — Jump in Pressure Coefficient Across a Flat-Plate Oscillating in Harmonic Pitch — 
M - 0.90. u> = 0.06 
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Figure 9. — Jump in Pressure Coefficient Across a Flat-Plate with Harmonically Oscillating 
Quarter-Chord Control Surface — M = 0.9, co = 0.06 
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Figure 10. — Jump in Pressure Coefficient Across a Fiat-Piate with Harmonicaliy 
Oscitiating Quarter-Chord Controi Surface — M =0.9, oj = 0. 120 
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Figure 1 la.— Jump in Pressure Coefficit 


Harmonic Pitch — M = 0 
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Figure 13. — Steady-State Pressure Coefficient for Airfoil Section 





Figure 15a — Jump in Pressure Coefficient Across an Airfoil with Harmonically Oscillating 
Quarter-Chord Control Surface — M = 0.875, cj = 0.05 
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Figure 15b.— Jump in Pressure Coefficient Across an Airfoii with Harmonically Oscillating 
Quarter-Chord Control Surface— M = 0.875, co = 0.06 



Figure 16a — Jump in Pressure Coefficient Across an Airfoil with Harmonically Oscillating 
Quarter-Chord Control Surface — M = 0.9, co = 0.06 
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Figure 16b — Jump in Pressure Coefficient Across an Airfoil with Harmonically Oscillating 
Quarter-Chord Control Surface — M = 0.9, ui = 0.06 






Figure 17b. — Jump in Pressure Coefficient Across an Airfoil with Harmonically 
Oscillating Quarter-Chord Control Surface — M = 0.9, oj = 0.06 


tion. is about the same as that shown by Ehlers (ref. 1 ) for lower Mach numbers. Thus the 
essential features of the experimental data are reflected in the calculations while the ampli- 
tudes from the two sets of data do not matcli well. 

Finally, figure 1 7 compares results from Freon and air calculations for the control surface 
case at M = 0.9 and co = 0.06. The Freon calculations for both the pitch and control surface 
motions reflect the slight forward shift of the shock shown in the steady results presented 
in figure 13. Except in the regions of the shock and hingeline, the Freon results also exhibit 
the slight increase in magnitude over the air results shown in the steady results. Surprisingly, 
the magnitude of the pressure coefficients in the neighborhood of the shock and hingeline 
appear less for Freon than for air. 

Generally, then, the results of this section are consistent with the results presented by Ehlers 
in reference 1 . The results for the flat-plate configuration correlate well with corresponding 
results from linear theory for both harmonic pitch and quarter-chord control surface motion. 
In all cases, the linear theory calculations provide slightly larger amplitudes than the finite 
difference calculations. For the airfoil, the correlation between the finite difference theory 
and the experimental results appeared to improve with Mach number. However, the calcula- 
tions, as in reference 1 , continue to provide pressure coefficient magnitudes larger than 
measured values except over the aft portions of the airfoil for the imaginary part. The 
correlation between theory and experiment remains inconclusive, and thediscrepencics may 
still be attributed to unknown problems associated with either the theory or the experiment. 
Finally, the set of two-dimensional examples presented here are limited in terms of fre- 
quency range due to the phenomenon of frequency limitation discussed at the beginning of 
this section, and also in section 4.2 and appendix D. 
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7.0 AN INITIAL APPLICATION TO A RECTANGULAR WING 


7.1 INTRODUCTION 

A pilot computer program has been developed for computing the unsteady transonic aero- 
dynamic flow over a three-dimensional rectangular wing. The program is based on the 
finite difference equations derived by Ehlers in reference 1 and represents a direct exten- 
sion of the two-dimensional program discussed both in that report and in the preceding 
sections of this report. 

A picture showing the geometrical setup is presented in figure 1 8. A rectangular, untapered 
wing is shown in the z = 0 plane, with leading edge at x = -1 and trailing edge at x = +1 . The 
wing tip is^t yo = y$> and a partial span control surface is included with a hingeline at 
X = Xa, an inboard side edge at y© = ya, and an outboard side edge at y© = yg- The steady- 
state velocity potential distribution, , is calculated for a wing of the same planform but 
finite thickness and evaluated at the finite difference points that span the mesh space. The 
flow is assumed symmetric with respect to the x-z plane at y© = 0, and thus the solution is 
carried out only over half the wing. The program is arranged so that mesh points lie in the 
plane at y = 0 and the boundary conditions of 

d\p 

— L = 0 

9y 

y =0 


results in (x,y,z) = (x,-y,z). This condition is readily included in the finite difference 

formulation. 

The equations used in the pilot program are given in appendix A. For the most part, the 
equations are directly from reference 1 . However, the second term in the expression for the 
velocity potential includes an integral in the flow direction with an upper limit of infinity. 
This term presents special problems in its evaluation. Morino et al.(refs. 1 2 and 1 3) present 
two evaluation procedures; a third is proposed here in hopes of obtaining improved 
efficiency. The formulation is presented in detail as part of appendix B. Briefly, here, the 
wake integral of the form 


OO 



^t(yi) 


d\p 


3x. 


f 


is converted to a form that may be evaluated using Laguerre integration. As in the two- 
dimensional case, advantage is taken of the pressure function 




3x 
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to evaluate far-field boundary conditions on the upstream and downstream mesh boundaries. 
Further, this expression is used on the upper, lower, and side boundaries to determine the 
velocity potential distribution on the boundary once a velocity potential has been deter- 
mined at one point on each flow-wise set of finite difference points. Again, for conven- 
ience, the velocity potential is evaluated on the boundaries only at x -values equal to the 
trailing edge (i.e., at x = +1). The .pressure function in finite difference form then permits 
the calculation of ip, for different! values of x for the constant values of z and y. 

Ehlers in reference 1 presents an alternative derivation for swept (but untapered) wing plan- 
forms. However, W. Schmidt (ref. 14) indicates that considerable success in the steady-flow 
problem has been obtained using a straightforward rectangular finite difference grid not 
necessarily aligned with the wing planform. This particular point would have to be investi- 
gated with respect to the unsteady transonic problem, but there is reason to expect the 
pilot program to be applicable with relatively little modification to planforms other than the 
unswept, rectangular planform for which it was developed. 

The major obstacle to performing three-dimensional solutions using finite difference 
techniques is the size of the complex velocity potential matrix, which must be stored 
between iterations. The finite difference mesh for practical cases is estimated to be on the 
order of 45 x 30 x 20, or some 81 000 words of core storage for the and matrices 
alone if the problem is to be stored in the machine all at once. The alternative is to store 
these matrices on tape and to bring in three (or more) planes at once. The current pilot 
program is capable of both modes of operation. 

Column relaxation has been selected as the solution procedure because (1) column relaxa- 
tion has proved most reliable for combinations of Mach number and frequency where the 
convergence is marginal, and (2) x-z planes require only three planes for the relaxation cal- 
culations of each plane. Note that x-y planes have this same feature while y-z planes 
require four planes because of the backward differencing of the x derivatives at super- 
sonic points. 

For our CDC 6600, a mesh of 25 x 1 9 x 20 is about the maximum size for an incore solu- 
tion. The 25 X 20 grid in the x-z planes appears small for practical problems, based on our 
experience with two-dimensional analyses. The problem size may be increased by going to 
an out-of-core program. Both the incore and the out-of-core versions may be run online. 

7.2 RECTANGULAR WING EXAMPLES 

The three-dimensional program was used to calculate the pressures over an aspect ratio 5 
rectangular wing undergoing harmonic pitch. Calculations were performed for both a flat- 
plate and a NACA 64A006 profile, with a Mach number of 0.875 and a reduced frequency 
of 0.06. A mesh of 44 points in the x -direction, 1 6 points in the y -direction, and 26 points 
in the z -direction was used. The extent of the finite difference solution area in physical co- 
ordinates was X =+3.80 (the chord ran from -1 .0 to +1 .0), yo=+10.5 (the wing tips were at 
+5.0) ,and Zq=+10.0. The computing time for these calculations ran on the order of 7 to 8 
sec of CPU time for each iteration and 8 to 9 sec for each far-field update. This was for a 
CDC 6600 computer using the KRONOS 2.1 operating system. The number of iterations 
required for convergence (in this case ERROR < 10“'*) was on the order of 1 80 when start- 
ing with zeros for the initial velocity potential field. 
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The flat-plate results are presented in figure 19. The chordwise distribution of the jump 
in pressure coefficient is prtesented for five spanwise locations. In each case, the finite 
difference results are matched wi^h results from the NASA subsonic program (refs. 4 and 
5). Correlation between the results from the two methods is good, and corresponds to that 
experienced with the two-dimensional examples of section 6.0 and reference 1 . 

The steady-state pressure distribution for a NACA 64006 profile configuration is shown in 
figure 20. It was obtained by using a program developed at NASA-Ames by Ballhaus and 
Bailey (ref. 1 5) with the mesh arrangement modified in the manner of Schmidt, Rohlfs, 
and Vanino (ref. 14). The corresponding jump in the unsteady pressure coefficient is 
presented in figure 2 1 for five spanwise stations. 

The results from a two-dimensional calculation are included in figure 2 1 a. Comparison of the 
two- and three-dimensional results reflect the anticipated softening of shock effects from the 
tluree-dimensional representation. In particular, the imaginary part of the pressure coefficient 
shows no apparent shock influences, and resembles what would be expected for the elliptic 
problem. Compare, for example, the linear three-dimensional solution shown in figure 49a 
with the linear two-dimensional solution shown in figure 4. Here the ACp changes sign at a 
significantly more forward chordwise location for the three-dimensional example than for 
the two-dimensional example. Three-dimensional experimental data are needed for further 
confirmation of the analytical program. 
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Figure 19d. — Pressure Coefficient Distributions for an Aspect Ratio 5, Fiat Plate in Pitch 
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Figure 19e. — Pressure Coefficient Distributions for an Aspect Ratio 5, Fiat Plate in Pitch 




Figure 20. — Steady-State Pressure Coefficient Distributions for an Aspect Ratio 5, Rectangular Wing 
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Figure 21b. — Pressure Coefficient Distributions for an Aspect Ratio 5, Rectangular Wing in Pitch 
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Figure 2 Id. — Pressure Coefficient Distributions for an Aspect Ratio 5, Rectangular Wing in Pitch 
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8.0 CONCLUSIONS 


Following the work of Ehlers (ref. 1), this report has further explored the characteristics of a 
finite difference solution to the unsteady transonic flow problem. Relatively rapid solution 
convergence rates were obtained using selected values of the overrelaxation factor together 
with row line relaxation. Additional examples were calculated for Mach numbers from 0.8 to 
0.9. Correlation of flat-plate results with corresponding results from linear subsonic programs 
was good for both pitching and control surface motions. Correlations between analyses for a 
NACA 64A006 airfoil and available experimental data were about the same at the higher 
Mach numbers, as reported in the preceding report (ref. 1) for the lower Mach numbers. Also, 
a pilot program for three-dimensional flow was developed and applied to a rectangular wing. 

Of significant concern was the encountering of an upper limit on reduced frequency, 
depending on Mach number and mesh region dimensions, above which the relaxation 
procedure will not converge. The relaxation procedure has been analyzed from several points 
of view, providing an explanation for the lack of convergence under certain circumstances. 


Boeing Commercial Airplane Company 
P.O. Box 3707 

Seattle, Washington 98124, July 1975 



APPENDIX A 


EQUATIONS FOR THREE-DIMENSIONAL FLOW 


The flnite difference equations for three-dimensional unsteady transonic flow are given by 
Ehlers in reference 1. They are listed here for convenience. At interior points where the 
steady flow is subsonic, we have 




= -E. 


’"Mjk+1 “> "^‘i+ijk ■ i-ijk * “Vj ■ S '^‘ii+ik 

and at points where the steady flow is supersonic, 

%'^‘ijk-i '(^yj '’yj "^^^k 


‘’^k ‘^‘ijk+1 ^ "^‘i-iik ‘ i-2jk ■ “yj "^‘u-ik * ^yj '^•u+ik 


where 


■% “ (Zk+i - ^k-i)U - ^k-i) “’'i ° 'yj+1 ■ yH><yj-yj-i> 


b,. =■ 


1 


b.. = 


1 


2k (zk+1 - 2k-iH2k+i -Zk) yj <yj+i • yj.|)(yj+i -yj) 


(Al) 


(A2) 


(A3) 


E, = 

Ej = di«i+>/2jk •'wd,./€ 

E 3 = Ci_| Uj .i/jjk-iwcjj/e 
E 4 “ di-l“i.3/2jk-iwd,i/€ 


“ (Xi+i • Xi.i) (Xj+^ - XjF = <’‘i'*i-l>‘^i 

“i+>Ajk = K - - ^oy,,)/(xH-l - Xi) 

“i-‘/^k” * ^®i-ljk^^*‘*^‘‘^^ 
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X: - X 


1. .. H-^i-2 


Equations (Al) and (A 2 ) are equations (24) and (27) from reference 1 . 

The boundary conditions on the upper and lower wing surfaces lead to the following 
equations for subsonic flow at finite difference points immediately below the wing, k = 




ijk. 


m 


= '^M+ljkn, -^.^.i-ljk^-^yj'^iij-lk^-^j'^.ij+lk^-hi S^Pij 


(L) 


(A4) 


and points immediately above the wing, k = kj^ + 1 

- (ay. + by. + + E, + E, - Qijkn^+l/^) 

"^li+ljkn^+l '^2^ii-ljkn,+l ‘Vj '^ly-.lkn^+l "'’yj y+lkn^+l %j^+l ^ 

where 

Ff) = ffW (Xj,yj) + l«fP (Xj.yj) 


(A5) 


(A«) 


The (L) and (U) refer to upper and lower wing surfaces, respectively. The equations simUar 
to (A4) and (A5) at supersonic points can be written down analogously. 

The total harmonic deflection of the wing is written as 

Zo = Sf(x,y,t) = 6|fo(x,y) + fj (x,y) e*^^^ I (^7) 

The steady velocity potential, <fio is calculated from the steady deflection shape, fo , while 
the unsteady potential ip ^ , is. calculated from the harmonic mode shape, fi (x,y). 

Over the wake, the conditions that the trailing vortex sheet supports no pressure, 


3x, 


ituAv^j 


= 0 


(A8) 
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results in a term being added to the right>hand side of equations (Al) and (A2). For finite 
difference points just below the wing plane (k = the additional term is 




(A9) 


and for points just above the wing plane (k = k„, + 1 ) the term is 


where 




Affi . .. = AsP . . 

' « * 1, + Ij 


,-MxpXi +i) 


(AlO) 


(All) 


and is the jump in velocity potential at the first point aft of the wing trailing edge at 

station j determined so as to satisfy the Kutta condition on t)ie trailing edge. The additions 
of equations (A9) and (AlO) implicitly satisfied the condition. The normal velocity is 
continuous across the wake. 


The finite difference equation for the jump in across the wing to the second order in 
mesh size is 


A<p^ 






where 


= 1 

^sl 4s, (s, + 1) 

h(2s, + 1) 
^sl ■ 4(s, + 1) 


(dsl Ff^-ds2Fj^^) 
1 

*^82 4sj(s, + 1) 

_ h(2Sj + 1) 

^s2 ■ 4(s, + 1) 

+1 


(AI2) 


(A 13) 


(A 14) 


Two integral relations are used to satisfy the far-fieid boundary conditions on the outer 
boundaries of finite difference mesh. The first for the velocity potential is 
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(A 15) 


■Vt 4(y;> 


-Yt Xt(y/) 

*^f 1 '<'xr‘"<^-'>'- '^^“x.'x.'l'"'' 

V 

and the second for the pressure function is 


P(x,,y,,z,) = + iw^j 


+yt /t(y/) 

i-j j [A,., X,;-XA^,^_,]dx,'dy.’ 

-Yt xj(y/) 


/yt 

J e*^*t^yi)A^,^(y,')^z, (My/)'’4>yi 'y/’^i^^^y/ 


-Yt 


■^45Firr V 

Jy ^ ^ ^ ^ 

t 2 2 2 

Defining x = x, - Xj', y = y, - y,', z = Zj - z/ and R =\“xj +"y +T 

_ iX, (MX - R) 

\Kx.y,z) = 5 -t: 


^z/ =-|-0/R + i^,)^ 


* = [iX,M -^(1/R + iX,)]^ 

X = +iw^ 

X = [iX,M-^(l/R + iX,) + iw] 

Xj.> = -^|[^ - iX, M - iwj (iX, + 1/R) - X" x/r| ip 


(A 1 6) 


(A 17) 


(A 18) 
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Equations (A 15) and (A 17) have been simplified for purposes of the pilot program. First,, 
as noted in the two-dimensional derivation 



(A19) 


and thus the second integral in the first term of both (A 15) and (A17) are zero. Second, the 
third term, which is the volume integral and has not been of significance in the two- 
dimensional problem, has been dropped. Third, since we are interested in the far-field, we 
approximate Xj - x/ with Xi and yi - y/ with yj so that the terms of and x may be moved 
outside the integral sign. The evaluation of the wake integral in equation (A 15) is discussed 
in detail in the next section. The equation for the velocity potential on the far-fleld (A 15) 
for Xi = 1.0 is 

+yt ;^t(yi'> 


/■'i 

I dx,'dy/ 

*yt ^«(y/) 


(A20) 


Where 


1 

I Vyi>z,;y;)dy; 

-yt 

and is defined in the next section. The pressure function (A9) becomes 

^yt *t(y/) 

P(x,,y,,z,) = I I A-Pjdx/dy; 

.r 


-yt xg(y/) 


Pt 


j e A^,^(y/)dy,' 


-yt 


(A21) 


(A22) 


Equation (A20) is used to evaluate the velocity potential along the line resulting from the 
intersection of the y-z plane through the trailing edge of the wing and the x-y and x-z planes 
bounding the finite difference volume. Equation (A 16) is then integrated by the trapezoidal 
rule to determine values ahead of this line and behind this line on the upper, lower, and side 
boundaries. For example, on the lower boundary when k = 1 , for points ahead of the trailing 
edge(x< 1.0ori<ii) 


’’■Mil " 
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and thb equation for points downstream of the trailing edge (ii > i) is 

'^‘Ul " + [P(Xi.yj.Zj) +P(Xi.i,yj,z.) e (^24) 

The application of equation (A22) to the upstream and downstream boundaries results in 
the following equations: on the upstream boundary 


*^*ijk ^Ijk 

and for thfe downstreaih boundary 


where 






Pljk = ^.(*1 >yj*='k) 


. (A25) 


(A26) 


and 


Pi ik = P^’^i yj’^k) 

max-*^ ‘max’ J 


max’ 


5, = X2 - X] 


62 = X 


1 • ’^i - 1 

‘max ‘max * 


(A27) 


1 + ia>6, /2 
‘^kl " 1 - ico5j/2 

Ck2 = 6,/(l -ico6,/2) 

Ck3 = (1 -iw5j/2)/(l +icj62/2) 

cjj 4 = 52/(1 +10)62/2) 

Equations (A25) and (A26) may be used to substitute for and (p, • in equations 
(Al)and(A2). ^ 
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APPENDIX B 


EVALUATION OF FAR-FIELD WAKE INTEGRAL 


The wake integral of concern, , is the second term of equation (A15), i.e., 


pi 


oo 

•I e (Xi -x,',y, -y,',z, -z,')dx/ 


xt(yi ) 

where the partial derivative is to be evaluated at z'l = 0 and ^ is defined in equation (A 18). 
The evaluation will be carried out for X| = 1.0 for a rectangular wing for which the trailing 
edge is Xj (y/) = 1.0. Equation (B1 ), after taking e into the Xj ' integral, becomes 

j Av?,t(y,')dy,'^ ® ■“^^’‘‘''*^'i^'U-Xi',yi-yi',z,-z, ')dx/ (B2) 

-Vt 1 

Let be the inner, x/, integral 

OO 

*w = y e-‘‘^(’‘i'-l)-||^,(l-x/,y,-y,',z,-z,)dx; (B3) 

1 

setting p = x/ - 1 and inserting the expression for ^ from equation (A 18) 

I - ? e-icop_3_ 

/ e dp 


withR, = Vp2 + Rg and R§ -(y,-y,')2 +(zi-z,')2 


Taking the outside and combining the p exponentials, we obtain 

OZi 


oo 

■if 


r M(w+X, M) 1 

■>mLmRi — »i 
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ior since 


/ «mA 

wM/ci-M*) / 

e'i'M (MRi +p) 

— 


MR, +p 

“ ^Ro /J ’ 

and as p-*> oo u -+ «> 

Further 

Vr+^ = -^[0^Ro" + M2(p?+Ro^) + 2MR,p+p^] ' 
which becomes, using = 1 - 

vTIV = ^ [ro^ +p2 +2MR.P +MV^] 


and 


so 




'/2 R, + Mp 


^Ro 


du L Tmo/r + il - ‘ = y/i±^ 

1F~ ^Ro + R. 0Ro Ri 


Xi/3Ro 


= -^ f^ULl 


du 


Next, let 


When 


M(w+X,M) 


Xi 


A 

*w dz 


U = 


oo 

r/ 


0 


u = 


p = 0 


(B6) 


(B7) 


(B8) 


(B9) 


(BIO) 


(BID 


(B12) 


(B13) 
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Now, let u = V + M//J so that 


oe , ^ 1 

3 r ■* M 

w = a^7 y -5 di^ 


0 Vl+(»'+M//J)* 

The singularities of the integrand arc where 

l+(i^ + M//3)2 = 0 


(B14) 


(B15) 


or 


which is in the left half of the complex plane. Thus applying Cauchy’s Theorem to the contour 
integral 


I 


w. 




i X,^Ro 
-i— («^ + M//J) 




c Vl+(»^ + M/0)2 
where the contour is shown in the following sketch 



Re(i^) 


and 


I = ***” 

^w r->oo 


(Bl6) 
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+ 

II 

o 

(B17) 

['3z,' *Cj ■ az,' *Ca] 

(B18) 

, = r e‘*^, 0 < 0 < 7 t/j , and 
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rg 


dg 
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*w “ * r-»-oo ai 7 


On C 3 , i» * ii? and di» = id »7 so that 



X.gR, 
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Vi+On+M/py 
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X.gR, 
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Integrating by parts with 
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(B27) 
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M 


d»7 


(B30) 


Note that while integrating by parts increases the order of the integrand, it reduces the number 
of numerical integrations from 4 to 2. 

•Anticipating the ultimate use of Gauss-Laguerre Quadrature, the following transformation is 
made; 

X,^Ro M _ 

•77 .or 77 = 


T =■ 


M 


Xj^Ro 


so that 


or 




(B31) 


w 
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where Zj ' is taken to be 0 in the evaluation of Rq 
L etting 


Kr)sV/J>+ (M-i^y 

the integral in (B32) has the form 

00 

J' f(T) e ■'*’ dr 


(B33) 


0 


This may be evaluated approximately as 


N 

i=l 

where w, and Ti are the weights and abcissas for Laguerre integration for a given N. 


(B34) 
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APPENDIX C 


TIME-LIKE CHARACTERISTICS OF THE RELAXATION TECHNIQUE 


Following Jameson (refs. 8 and 9), we study the convergence of the difference methods by 
considering the iteration numbers as a time-like variable. Thus the relation between con- 
secutively iterated values at an arbitrary point ij is defined in terms of a time derivative by 




(Cl) 


where the superscript denotes the number of the iterations. 

Cl COLUMN RELAXATION WITH INCREASING i 

The difference equation for column relaxation at subsonic points with i increasing is given for 
two dimensions by modifying equation ( A 1 ); namely 



where the superscript s denotes the value obtained from the solution of the equations for the 
ith column and n and n-1 denote updated and non-updated values, respectively. 


When we apply a relaxation parameter r, then we obtain a modified value of according to 
the formula 


Cs) 

Solving for - 






.(nrO 


Finally, eliminating y by means of equation (Cl) leads to 


(C3) 


= v»(^> At • (C4) 

To obtain the differential equation corresponding to the difference equation (C2), we introduce 
equations (C4) and (Cl) into equation(C2) and take the limit as Ax and Ay go to zero. We 
have 
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- Xj.j ) (Xj - 

Xi.i) ■ 

' (Ax, 

1 + AX2>AX2 



(C5) 


(C6) 


Substituting equation (C6) into equation (C5), expanding u and (^, about the point ij, and 
then simplifying, yields 


(U ^ 

2io) ^ , . 

- g + 

V?iyy + qv’i -2u V?,xt 

JL 

/ At V 

P - <*•-!) , 


Ax, 

\Ax, +AX2/I 

Lr r 



At \ri/^w 

(^*0 + ntAr Av Av^ 

^ e * 

kAx, +Ax2 / LrxAxj/ 

r Ax2j'^‘t + 0 (At,Ax,Ay) 


(C7) 


Since the real part of the coefficient of (^, ^ must be negative to produce damping, then we 
must have 

J_ ( r - 1 ) I 
r ' r VA \2 ! 
or 

^ 1 + Ax, \\ 

Axj/Ax, 

In practice, values of r in excess of this limitation have worked sir essfully. For example, 
with a spacing factor of 1.4, 1 + Ax, /Ax, = I 71 '’ ! ’ ■ ii, > w 'oilu i ihe sweep 

direction is in the direction of increasing or uecie..Mi:g .n, m. m. . , , ii. .e su^cesslully used 

values of r up to 1.9. At elliptic points, we overrelax and hence r is restricted to values 
between 1 and 2. 


It was suggested that complex values of the overrelaxation parameter might be useful in 
speeding convergence and a few values were tried without success. We will investigate the 
effect of a complex value of r on the real part as the coefficient of the >p^ ^ which provides 
damping. In equation (C7), the coefficient of the first term becomes 
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-2vUX2 



= -2ua 

2 




where re*'^ is the complex relaxation factor. The real part of the coefficient provides the 
damping and must be negative. 

(1 +8) COS 7 - r ^ Q 


We see that the imaginary component always reduces the amount of damping. Since generally 

u/AX| > to/e 

the first ^ provides most of the damping. With a complex relaxation parameter, the term 
proportional to 2ico/e also may contribute to the damping. The coefficient from this term 
becomes 

The real part is seen to be 

IcjOOCi (l-6|^)sin7 

e5r 


This term provides damping when 

5 < 1 for which 7 < 0 


or when 


5 > 1 for which 7 must be > 0 


Since 7 is a constant and 6 takes on values both greater and less than 1 , the effect of the 
imaginary component of the overrelaxation factor is to reduce damping; and hence no 
improvement in convergence can be expected from the use of complex values. The few runs 
made with complex relaxation factors support this conclusion. 


For convenience, equation (C7) may be written 


(u<P, x>x - *^ 1 X yy 'P, xt 


- 2 uo £2 'Pif + licoa^^Pn = 0 


(C 8 ) 
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where | 


a, = 


At 

Ax, + Ax 2 




Equation (C8) is hyperbolic in time and this is readily seen by eliminating the y?, 
derivative by a new time variable 


T = t + a, X 


This leads to 


2ico 


<“^x>x * '^•yy '^•tt 

-2(ua2 - icjttj)^,^ = 0 


(C9) 


Since u is positive, the coefficients of ip, and (p, ^ have the correct signs for a damped wave 
equation. 

We now consider hyperbolic points. When the local flow is supersonic, backward differences 
for the X derivatives are used to eliminate downstream influences. The difference equation 
from equation (A2) is given by 

2<=H “M (»!? -■'ti-lj) 

■ <■. M f (CI0> 

Using equation (A3) for the definition of the coefficients and eliminating tp\ by equation 
(C4) yields the following differential equation; 


(U<P, X + *^1 


yy 


+ q'/’i 


“Ax, (Ax, +Ax2> ■ e (ax, "^Ax, +AX2) ( r 


(Cll) 


where 


' Ax, = Xj -Xj.j , Ax 2 = Xj.j - Xj _2 


This equation is of parabolic type. Since u < 0, we must have r < 1 to have damping. Thus, at 
hyperbolic points, we use an underrelaxation parameter. 
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SimUar results are obtained for column relaxation with decreasing i (upstream).'In place of -, 
equation (C8), we obtain 

(u«p, x>x - X yy xt 


where 


-2 (utt4 - icotts ) V’l t = 0 

^ AXi [_r r AX| J 

\ 


The transformation t = t - ajx yields a hyperbolic equation of the same form as (C9). For 
supersonic points sweeping in the direction of decreasing i, the superscript (n) in equation 
(CIO) is replaced by (n - 1). The time-dependent differential equation then in place of 
equation (C 1 1 ) is: 

(uv^i — 'Pi ^ + 'Pi yy + q >Pi 

At \ ^ . 2 \ o 3 1 . 1 \ _ „ 

\AXi (Axi -HAxj)/ r e \Axi Axj +Ax 2 /'^‘t 

Decreasing r serves to increase the damping in supersonic regions. . : • 

C.2 ROW RELAXATION 

For the smaller frequencies, the use of row relaxation improves the rate of convergence. 

This may be attributed to the coarse mesh in the y direction. For row relaxation in the 
increasing j (or y) direction, equation (C2) for elliptic points takes the form 

Uj+ij (»^f iiij -V’i"]) -2diUy(.pf®^^ - 

(C14) 

+2aj(4'y-l - V>f^)-2bj(.pl®j - v?!") Vl) + Qy = 0 

Substituting equations (Cl) and (C4) into equation (Cl 4) and taking the limit as Ay, Ax ^ 0 
yields 

(u«Pix)x ■2-^v’ix ■'■'^lyy ''■‘IsPi 


(Cl 4) 
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where ^ and bj are defined in equation (A3) and Ayj = yj^.j - yj. j. Since the coefficient of 
^ must be negative,. this would appear to limit the range of overrelaxation factor to 


ai + bi Vi+l-Vi-l 

r<-J i= 1 + Li < -7 

aj yj-yj.j - 


The transformation 




, r 

yields the equation 

icjj r a; + b; 1 

(uV’ix)x -2— «^ix + ^lyy + 

-(At/Ayj)^ = 0 

which is hyperbolic in t. 

For supersonic points, we write in place of equation (CIO) for increasing) 
2ch Uj| ( 4 *^ - ) ■ 2di.j Ui.jj 


(Cl 6) 


(Cl 7) 


= 0 


Substituting equations (Cl) and (C4) into equation (Cl 7) yields 

(uv?ix>x ■'■'^•yy 

Since u < 0, this equation is hyperbolic with respect to y since the coefficients of and 

ffi. terms will both be negative and opposite to that of v?, . To render the x variable 

X X y y 

time-like, we write the equation in the form 


(u x) 


2ico 


x'x 


^Ix *^iyy ■*■^^1 


[v’lyt ^2 •Pit] “ ^ 


(CIS) 
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and deteimine^t and to eliminate the <fii yj, ^ ^ tenns and to obtain a coeflicirat 
of ^ witii the correct signl Choosing r = t I y t o(| X we And that- 






= 


(u^ - 2iw/e) 
2u 


Substituting «i > fit relation for r, we obtain the following differential equation 

equation for ^ 


(“'^»x)x - ^^»x + '^«yy -(^j) 0^'" “ 0 

In Older for the coefficient of the term to be positive, we must have 

>-u 

0, >\/iur 


(C20) 


In the light of equations (Cl 7) through (C20), the differential equation to be differenced U 
2i<u 


(u^ix^x * £ X ^ ^^*yy 
LX 


At , r (u„-2iw/c) 

-2*2^7 c ViuT [<Axt + 2 u '^‘t] 


(C2I) 


-2aj(-^^) At^it = 0 

where c is a constant > 1 . The difference relations for ip, ^ ^ and ^ are seen to be 

■ ^5 ?” * /.h/) 

The addition of the time-like differences caused the relaxation method for mixed flows to 
converge, whereas the difference equation (Cl 7) caused the iteration to diverge. 


(C22) 
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APPENDIX D 


I 


MATRIX ANALYSIS OF RELAXATION TECHNIQUES 

As noted in the text, in the course of the solution of the unsteady transonic problem using 
line relaxation (block successive overrelaxation, or BSOR), it was found that convergence 
was obtained only for relatively low values of w, the reduced frequency. In this appendix, a 
derivation is given that provides the theoretical basis for this empirical observation. By 
simplification of the basic boundary value problem and use of a theorem on the convergence 
of relaxation procedures, the conditions are determined for which the system matrix ceases 
to be positive definite, and thus for which a necessary condition for convergence is violated. 
The condition is defined in terms of the dimensions of the solution region, the transonic 
similarity parameter K, the Mach number, and the reduced frequency. It is convenient to 
think in terms of a critical reduced frequency, below which the solution procedure 
converges and above which it diverges. As shown in section 2.2 the co^p as derived here, 
corresponds very closely with those values encountered during the running of sample problems. 

The two- and three-dimensional problems will be developed together. The mesh regions for a 
two-dimensional problem is a rectangle with sides of lenths a and c. The region for the three- 
dimensional problem is a rectangular parallelepiped with the additional spanwise dimension 
of b. In order to carry out the convergence analysis and obtain an analytic formula for cu^p 
we make the following simplifying assumptions: 

• The steady-state perturbation potential, , is constant; 

• The far-field boundary values are constant; 

• The airfoil and wake boundary conditions are omitted; 

• The mesh spacings resulting from the discretization are uniform in each of the coordinate 
directions; 

• The relaxation factor is the same for all points and for each iteration. 

These assumptions, particularly the third, are rather strong; nevertheless, analysis of the 
problem resulting from these simplifications is found to yield results in generally good 
qualitative and quantitative agreement with experimental observations, indicating that the 
essential frequency-limiting characteristics of the problem have been retained. 

In mathematical terms, the effect of the assumptions is to reduce the problem to the solution 
of 

Kv^ixx +-/’.yy - ^^ix = 0 

over a rectangle (or the analogous differential equation and region in three dimensions) with 
Dirichlet boundary conditions, and it is this problem to which we apply the convergence 
analysis. 
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Following the sections giving mathematical background and details of the region discretization, 
we give the analyses for both the exact problem above and a modified problem resulting from 
transformation of (Dl). The rather complicated formula for cogj. from the exact problem is 
seen to reduce to the much simpler formula for from the modified problem as hj^ 0. 

D l MATHEMATICAL BACKGROUND 

In this section we give a brief statement of various mathematical facts that will be required in 
the course of our analysis. 

D.1.1 MATRIX FORMULATION OF BSOR 

Our discussion is a specialization of that given in Varga (ref. 16), and the convergence theorem 
we state is an immediate consequence of the one given there. 


/ 


Let an n X n Hermitian matrix A be partitioned for some N, into the form 



(D2) 


where the diagonal blocks Ajj, i = 1, . . ., N are square and nonvoid. Then defining the n x n 
matrices D and F by 



and with F* as the complex conjugate of F, so that A = D -F -F*, then BSOR applied to 
A^ = R, where ^ and R are n dimensional vectors, may be written i 

(D - rF)i^"»+^> = {rF * + (1 -r)D} + rR (D4) 

m = 0, 1, . . ., where r is the overrelaxation factor. The following theorem gives necessary and 
sufficient conditions for convergence of the above iteration. 


Theorem 7.— If D is positive definite and 0 < r < 2, then the given iteration is convergent for 
all starting vectors if, and only if, A is positive definite. 

Comments: 1 ) An implication of the theorem is that under the given conditions the eigen- 
values of the iteration matrix (D - rF)- { rF • + ( I - r)D) 
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will all be less than I in modulus, a well-known necessary and sufficient condition for 
convergence; and 2) the statement ''all , . . does not limit the applicability of the 
theorem. In any practical problem a stationary iteration such as (D4) must converge for 
every choice of starting vector to be of use. 

D. 1 .2 EIGENVALUES OF A CERTAIN TRIDIAGONAL MATRIX 


The eigenvalues of the tridiagbnal matrix 


a b 

c a 





are given in Bellman (ref. 1 7) by 

Xj = a - 2 Vbc ®os~^, ft — 


n. 


(order n) 


(D5) 


D. 1 .3 KRONECKER (TENSOR) PRODUCT 


From Bellman (ref. 1 7) we have the following: 


Definition.— IM A be an N| dimensional matrix and B an N] dimensional matrix. The 
Kronecker product of A ahd B, written A® B,.is the N| *N 2 dimensional matrix given by 


*21 B *2 2® 


*iN,B 

*2N,B 


and 


[aN,«® “N,2® 


"N.N, 


B 


Theorem 2.-Let Xg, £ - 1 , . . ., N| be the eigenvalues of an order Ni matrix A, and let 
1, . . ., Nj be the eigenvalues of an order N 2 matrix B. Tlien the eigenvalues of the 
N| *N 2 order matrix 

ONa* *N,) 


where I denotes the identity matrix of the given order, are given by 

Xg + gij for £» 1, . . ., N| and 1, . . ., N 2 . 

Comment: If B is tridiagonal, then the matrix (I}y|^ ^ A) + (B ^ is easily seen to be 
block tridiagonal,, which is the case we will encounter. 
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D.2 REGION DISCRETIZATION AND MESH POINT ORDERINQ 
D.2.1 TWO DIMENSIONS 


In two dimensions the region is a rectangle in the x-y plane with sides of lengths a and c in 
the X and y directions, respectively. Taking the lower left comerpoint as (x©, y©) and 
chooang N^ and Ny as the number of mesh increments in the indicated directions, the 
totality of mesh points is given by (xg, yj), where 


Xg = x^+£hjj, £ = 0,1, . . „ Njj, 

yj ~ Yo Jhy, j = 0, 1 , . . Ny, 


(D6) 


with hj^ = a/Nj^ and hy = c/Ny. 

Of these, points with £ = 0 or N^, or j = 0 or Ny, are boundary points; all other are interior 
points. Thus the 

Ns = (Nx-l)(Ny-l) 

interior points are given by equation (D6), but with the ranges of £ and j being 1 to N^ - 1 
and 1 to Ny - 1 , respectively. 

Now to obtain a matrix formulation of the discretized problem, the interior points must be 
ordered as a one-dimensional array. To be definite, we choose column ordering, which is 
consistent with column relaxation of the resulting system. A similar analysis may be done for 
row ordering, which is consistent with row relaxation and leads to precisely the same formulas 
for tOgj.. The column ordering is specified by having the j index vary most rapidly and £ 
index least rapidly. That is, by the mapping: 

(X£, yj) pg, where S = j + (£ - 1 )(Ny - 1 ) (D7) 



92 




D.2.2 THREE DIMENSIONS 

In three dimensions the region is a rectangular parallelepiped with sides of lengths a, b, and c 
in the x, y, and z directions, respectively. The discretization in the x and z directions is the 
same as that given by equation (D6). We will treat the y (spanwise) direction slightly 
differently since we have a symmetry condition in the x-z plane. In particular, taking yi = 0 
and yj»j = b/2, the totality of mesh points is given by (xj, yj, zj^), 

xj = Xo + Ehjj, C = 0,1 Njj, 

where 



II 

j= 0,1, . . ., Ny, 

(D8) 

and 

Zy. = Zo +khjj: 

, k = 0,l N^, 


with 

hx = a/Nx, 

hy=b/2(Ny-l), 


and 

hz 

= c/N,. 



Of these, points with C = 0 or Nj^, j = 0 or Ny, or k = 0 or N^ are boundary points; all others 
are interior points. Thus the interior points are given by equation (D8) but with the ranges of 
£, j, and k being 1 to N^-1, 1 to Ny-1, and 1 to N^-l, respectively. 

The three-dimensional mesh point ordering is a direct extension of that used for two 
dimensions. That is, the interior mesh points are sequenced such that the k index varies most 
rapidly and j index least rapidly. Thus the ordering is given by the mapping 

where (xg, yj, zj^) ->-PS, 

S = k + (g- DCN^-l) + (j-l)(Nx-l)(Nz-l), 

forC= 1 Nx-lJ= 1 Ny-1, 

and k = 1 , . . ., N^-l. 

D.3 DERIVATION AND ANALYSIS OF THE SYSTEM MATRIX 

D.3. 1 MODIFIED PROBLEM-HELMHOLTZ EQUATION 

Since it will be most economical to derive the system matrix for the modified and exact 
problems at the same time, we first define what we call the modified problem. This consists 
of the given region and the differential equation obtained from equation (Dl) or its three- 
dimensional analog by making the transformation 

(Pj t// e^^x/eK 
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I III mill ■mil I 


iiiiiii iiii iiini II I II II II I 


I II I I I II 1 1 


When this is done we have the equation 


^ ^VV 






\p = 0 


(DIO) 


in two dimensions and the equation 




XX 






zz 




^ = 0 


(DU) 


in three dimensions. We observe that the first derivative term is absent in these equations— a 
consequence of, and the motivation for, the given transformation. Further, we observe that 
except for the K coefficient, which may easily be removed by scaling, (DIO) and (D1 1) are 
the two- and three-dimensional versions, respectively, of the well-known Helmholtz (or 
reduced wave) equation. Thus the modified problem is directly analogous to the mathematical 
formulation of the problem of a finite elastic membrane whose amplitude is assumed to vary 
sinusoidally in time. 


D.3.2 TWO-DIMENSIONAL ANALYSIS 


Note that both equations (Dl) and (DIO) have the form 

Kip^^+>fiyy-2iCi^^ + C2<p = 0 


(D12) 


where ~ (DO 

^2 

and Cl = 0, Cl = for equation (DIO) 

e(l-MO 

From this fact it follows that the system matrices from the discretizations of (Dl) and (DIO) 
have the same form, which we now derive. 

At each of the N^ interior points of the mesh given by equations (D6), we discretize (Dl 2) by 
making the substitutions 

*^'^xx^('^«-l,j +V’c+i,jyhx 

j-1 + n, j-u) /hy 


'^x“^('^)e-U, j j)/2hx> 


(D13) 


.where hj^ = VlC . 


This yields the system of difference equations 


- 82 ^ 8 - 1 , j -8i'^8,j-l + 6o'^«,j 


•8i'^K,j+l •fev’g+ij 


= 0 , 


(D14) 
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for K = 1 , . , Njj-1 , and j = 1 , . . Ny-1 , where 


go=2/h| +2/h* -Cj, g, = 1/hJ 

and g] = 1/hj^ + ici/h^. Here terms for which the first subscript becomes 0 or Nj^, or the 
second becomes 0 or Ny, are transferred to the right-hand side. 

Using the mapping given by equation (D7), equation (D14) may be written as 


S= 1,2, 


■82'^S-(Ny-l) 

-g j ^S+ 1 ■ 82 **^S+(Ny- 1 ) 0, 


(D15) 


N 3 , with terms being transferred to the right-hand side where appropriate. 


Now since each subscript is linear in S with a coefficient of 1, it is clear that (DIS) represents 
a five-diagonal system for the unknowns Letting‘s be a vector with these 

components, (D15) may be written in matrix notation as R, where‘s is the right-hand 
side vector and A is the Nig order system matrix given by 


where 


and 






A 




go 

-gi 


So -g, O 


Go - 



go -gi 
-gl go 


(order Ny-1), 


Gl =-gj *Njrl- 


(D16) 


(D17) 


Now we have the system matrix; it remains to apply theorem 1 to obtain an analytic form for 
First, we identify A as given by (D16) as a matrix of the form given by equation (D2) by 
rtiaking the identifications 


NfNx-l;Aii = Go,i= 1,...,N; 

^i+1, i “ ^*i, i-H “ Gl , i = 1 , . . N-1 ; and 
A^ = 0 otherwise 
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Next we check that D is positive definite. From the identification, we have that D as given 
in (D3) is 

D = diag (Go, Go, . . Go). , 

Thus the eigenvalues of D are the same as those of Go , which from (D17) and' using (D5) are 

go -2g, cos-j;|^,£ = l,. . ., Ny-1. 

The minimum eigenvalue occurs when £ = 1 , so that the condition on D is satisfied if, and 
only if, - 

go -2gi cosTjp>0 

I>y 


or, using the definitions of go and gi , if, and only if. 


2^2/, 7T 

-=— + ( 1 -COST?- I > Cj 

hi h= ( Ny/ 


(D18) 


where Cj =— for the exact problem and Ca = 




e(i-M 2 ) 


for the transformed problem. 


Condition (D18) will be satisfied, as will be apparent (eq. D19) since the frequency for which 
A ceases to be positive definite is substantially below that for which condition (D18) is 
violated. 


The eigenvalues of A are easily found through the use of the Kronecker product. Noting first 
that A may be written as 


where 


(a " *Nj^- 1 ® ® Ijvjy. j ^ 


O -gt 
-E2 O -g* 




- 82 ' o’ -g* 

-82 O 


(order 1), 


we have using equation (D5) and theorem 2 that the eigenvalues of A are given by 

_ £7T . / * jTT 

80 -2gt cos-TT- - 2 Vg 2 g 2 cos^ , 

I'ly 

where 

£=!,..., Ny- 1 and j = 1, . . ., N^- 1. 
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The minimum eigenvalue occurs when S = j = 1 so that A is positive definite if, and only if, 

go -2g, cos jj- - 2 y /%2 gj cos -jf" > 0, 

y ^ 

or using the definitions of go , gi and g 2 if, and only if 


2 


K 




(D19) 


Now we consider the two cases corresponding to the two different choices for C] and C 2 . 




Case /.—For the modified problem, Ci = 0 and C 2 = ' , 

e( 1 - ) 

and hence, by theorem 1 , BSOR converges if, and only if, 


so that A is positive definite. 






eK(l -M2) 


or if, and only if, w < Wpj. where 


0 ) 


M 




•/2 


(D20) 


Equation (D20) gives the precise value for in the transformed problem. A more 
transparent formula may be obtained by using the first two terms of each cosine series 
expansion and the relations hj^N^=a, hyNy = c so that to a very good approximation for 
reasonably large and Ny. 


CO 


cr 




*/2 


M 


b 


Kc- 


(D21) 


The dependence of co^j. only on the absolute dimensions of the region and not on the mesh 
size is notable, and has been previously observed in another context by Molberg and 
Reynolds (ref. 18). 


Note: In going from equation (D19) to equation (D20) we have used the relation 

l/eK = M2/(l -M2) (D22) 

2 

Case 2.— For the exact problem, and , so that, making these substitutions 

in equation (D19), we have that BSOR converges if, and only if, 



1 



(D23) 
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/ 


We note that this inequality is satisfied for cj = 0 and that as io increases, the left-hand side 
is strictly decreasing while the right-hand side is strictly increasing, so that there is a unique 
value at which the inequality just fails, namely 

An explicit formula for co^j. may be derived from (D23). When the > is replaced by an equal 
sign, (D23) may be transformed into a quadratic equation for which is easily solved. 


When this is done it is found that the solution may be written as 


where 


and 



(D24) 


with 


I 

2cos(7t/N^) / 

[i -2M' ^sin^^^^j+^Qi (1 -M^)M2 h| + [l -2M2^sin ^ 


We observe that as h^ ^ 0, Qj ^ 1 , so that ^ which is the same as equation 

(D20), the result for the transformed equation. 


D.3.3 THREE-DIMENSIONAL ANALYSIS 


The three-dimensional analysis is a straightforward extension of that already done for two 
dimensions. First we note that the three-dimensional analog of equation (Dl) is 

^2 




1 XX 




yy 




2ico 


zz 


■'Pi' 


. -A 

+— sPi - 0, 


so that both (Dll) and (D25) have the form 
where, as before 


'^yy "^zz + C2<fi = 0, 


(D25) 


(D26) 


Cl =0, C 2 = /e( 1 - ) 


for equation (Dll) and 


C] = oj/e, C2 = o)^/e 


for equation (D26). 
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At each of the Ng interior points given by equation (D8) we discretize (D26) by making the 
substitutions 

^yy-"('^Kj-l,k -2^fiJ,k +Vfij+i,k)/h2, 

'^zz‘*’('^«,j,k-l *2^fi,j,k +'^fij,k+l)/*^z’ (D27) 

"^x “*■('^8+1 J,k - '^8-1 J,ky2hx> and 

'^■^'^8,j,k. 


This yields the system of difference equations 

-«3V’«,j-l,k -feV?g.ij,k -8iV>fij,k-l +8o«^>fiJ,k 

-Ki'^gJ.k+l - 82 <P 8 +l,j,k ■fe'^BJ+l.k = 0 


for 8 = 1, . . N^- 1; j = 1, . . Ny - 1; and k = 1, .... N^- 1, where g 2 is as before, 
while now 





■+— -Cj , g, = — , and gj = 




(D28) 


Here terms for which the first subscript becomes 0 or Nx, or the second becomes Ny, or the third 
becomes 0 or N^, are transferred to the right-hand side. The ^gj.jk term, for which the second 
subscript becomes 0, is handled according to the symmetry condition: 

•^8,0,k = '^8,2,k 


Using the mapping given by equation (D9) equation (D28) may be written 


forS= 1 


•g3<^S-N2 ■*2'^S-N, *8i'^S-I ■•■So'Ps 



-gi'^S+1 *82'#’S+N, ■fo'^S+N2~®’ 


(D29) 


where N, =N 2 -l,andN 2 = (N^- IKN^- 1). 

Equation (D29) represents a seven-diagonal linear system, which may be written, as before, 
as A^= R, where R is the right-handi-side vector and A is the Ng order matrix given by 


G4 2G3 
G 3 G 4 G 3 



A = 



G 3 


G 4 G 3 
G 3 G 4 


(D30) 
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where and is the order matrix given by equations (D16) and (D17) with, 

however, the redefinition of go for the present three-dimensional case. 


We have the system matrix; now we apply theorem 1. With regard to the hypothesis, we 
observe that the modus operandi is stUl column relaxation, so that matrix D is still given by 

D = diag (Go , Go , . . Go ) 

and thus, following the previous argument, D is positive definite if, and only if. 


Til h* 




•y ..j, 

which will hold for the frequencies of interest. 


(D31) 


The eigenvalues of A may be found, as before, with the aid of the Kronecker product. We 
note that A may be written as 


where 


^ - (*Ny-l ® ® *Nj) 


O -2g3 

-S3 o -gj Q 



O 

-g3 


-g3 

O. 


(order Ny-1) 


Now the eigenvalues of G5 are given by 

-2g3COs^^, j= 1,2,..., Ny-1, 

while those of G4 are the same as those for the two-dimensional A, with go redefined. 

Thus, applying theorem 2, the eigenvalues of A are given by 

go -2g, cos^ -2 Vgig? cos||^ -2g3 cos , for 


S = 1 , . . ., N„*l ; j = 1 , . . ., Ny-1 ; and k = 1 , . . ., N--1 . 


The minimum eigenvalue is clearly that for which £ = j = k = 1. Thus using the definitions of 
Ko . 8i . 82 . and g3 , and applying theorem 1 , we have that BSOR converges if, and only if, 




Kh2 


(D32) 
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which is the three-dimensional extension of equation (D 1 9). 

Specialization of equation (D32) for the two cases follows directly as before. 

Case 1.— Fox the modified problem we have that BSOR converges if, and only if, < cUd- 

wliprp 

^ Lhx^ ' Nx / Kh* ' 2(Ny-l )' 


CO 


cr 


. 2 /, It 0 '^ 

H ( 1 - COS Tt— )1 

Khi ^ J] 


(D33) 


or, for N^, Ny, and N^, reasonably large, 


CO 


cr 


M La^ Kb* Kc*J 


»/2 


(D34) 


Case 2.— For the exact problem we find, proceeding as in the two-dimensional case, that o>d 
is again given by equation (D24) with Q 2 and Q 3 defined as following that equation, but 
with redefined as 


Qi = (l - cosi^^ +— ^ — (1 - cos \ • cos-^\ 

h* V Kh* V 2(Nv-l)y Kh* V V 


Kh* 


Kh| 


(D35) 


Thus, as before, the result for the exact problem reduces to that of the transformed problem 


as hx 0. 


D.3.4 SUMMARY 


In summary, the predicted values for resulting from the analysis of the simplified 
problem are given by equations (D21) (approximate) and (D24) (exact) for two dimensions, 
and by equations (D34) (approximate) and (D24) and (D35) (exact) for three dimensions. 
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APPENDIX E 


VON NEUMANN STABILITY ANALYSIS 


One way of investigating the stabflity of a differencing system is to study the propagation 
of errors through the mesh in the iteration process in a manner derived by Von Neuman (ref. 
19). Let E|^£ be the error in ifii ]^£ at the k, £ mesh point. Von Neuman assumed that these 
errors may be expressed in the form 

P“‘*^ax 9 ^*max 


for a uniform mesh. Here a and b are x and y dimensions of the mesh. Since there are the same 
number of equations as unknowns, the system of equation (El) is determinate. For simplicity, 
we shaU assume a flat plate with a uniform mesh. Then uy = K and the difference equation for 
the column relaxation sweeping in the direction of increasing i is 



(n) 

M-lj 




+ ip 


(n-1) 

•i+lj 
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/Ax^- 
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eAx 




+ 



y u-lj ^ y 


(E2) 


Now 


and we obtain 
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+ <fi 


,(n-l) 

i+lj 



i<o 

eAx 



- (1 -r) 
r 



<«)^[ 


^*ij+l ‘ij '^»ij-lj 


/rAy^ 


y+1 *ij 


!l'.V 


r e V r /^ly 




(E3) 


Since the difference equation is linear, we may consider a single component of equation 
(El) and find the solution of the difference equation in the form 
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e 


(E4) 


This has the initial condition in the form of equation (El). The iteration will converge if the 
magnitude of g is less than 1 . Substituting equation (E4) into equation (E3) yields, after 
dividing out the common factor, 

- ^e*^' -ge + 2g(cosflj - l)/r-Ay* (E5) 


CO^g 


(-!-^)(COS02 - 1) + 

^y 2 V r 2 er 


€ 


ef-y-o 


where 01 = 2itp— and 02 = 2n q— . Solving for g and simplifying, yields 
a c 


g=' 


A - c/r 


(E6) 


where 


A = a| + ibi 



2K 2(1 -COS02) 

Ax^ Ay^ ^ 


(E7) 


co^ 2K 

The quantity C is always positive if ~ < — • Since the magnitude of g must be less than 

one, we must have 


[ai +“ cj +b] < (a, --IrC) +b^ 

which simplifies to 

(2^)C (2a, -C)<0 

Since C> 0 this leads to 
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/KCOS0, <0 . ^ 

21 — +"r7rrsin0 

\Ax* eAx 


.V— 

/ Ax^ 


2(1 -cos^a) cj2 
Ay^ ® 


(E8) 


Considering the preceding relation as a restriction on frequency, we express the inequality in 
the form 


2 -L 


00^ + 


Ax 


r 1 - COS01 1 - cosOj”! 

sine j < 2Ke + K-Ay^ J 


(E9) 




or for M near the sonic value 


For p = q = 1 the limit Ax, Ay -*• 0 yields 

- <4K€ , 

V Kc^ 

M= \a’ KcV 

1 - 

since = Ke < 1. Note that this limitation depends upon the outer dimensions of the 

mesh region and is similar to the result from the matrix method. 

For now relaxation, with increasing j, we have 


O) 


(ElO) 


\^^‘i-lj ^>ij ^'i+ljj' V ‘i+lj ^'i-iy 

+ k":! ; -2^?.^. + , 1 /Ay^ = 0 
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(Ell) 


Eliminating if> 


,(s) 




(El 2) 
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Substituting equation (E4) into equation XE 12) and simplifying, yields 

(“"*.)[? ■ (t^)] 

[t -2(^) 

Solving for g yields 


-id 




A-j-C 


(E13) 


where 


A^e-'^'VAy 


.T = J«V 


^A = e VAy*;C = 


,_ 2 _ 

Ay* 


2K(l-cosd,) 


Ax? 


w; 

€ 




The condition then lgi<l becomes 


2cosdj/Ay* 



2K(l-cosd,) 

Ax* 


CO* _ 2<o 

c ” €*Ax 


sind, 


or 


2<o - 2K( 1 - cosd I ) 2(1 -cosdj) 

— + Sind, <■ - + —■ Ayi““ 

which is the identical condition found for Column relaxation. 


Since complex values of ORF were tried, it is worthwhile to find the effect of complex r on 
the Von Neumann test for convergence. Equation (E6) for g becomes 
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and the condition that lgl< 1 becomes 


+ + <(a-i^c)‘ +(b+-^)‘ 

which simplifies to 

(2a-C) <0 

For small frequencies, the same inequality (eq. E8) holds if 

r<2cos7 

Otherwise Ig l< 1 may hold for large values of frequency (for r > 2cos7). 

Oswatitsch and Singleton (ref. 20) replaced the actual time derivatives in the Euler 
equations of fluid flow with fictitious time derivatives to change the equations from hyper- 
bolic in time to parabolic in order to improve the convergence rate for calculating steady 
flows. In the same manner, an artificial time derivative was added to the unsteady differential 
equation for harmonic motion and an ADI (altemating-direction-implicit) scheme was tried 
to solve the resulting equation. For the flat plate, the equation is 


„ 2itJ , . w 

dt ixx e *x *xx e 1 




(E14) 


The mesh was swept through by row relaxation and alternately by column relaxation. 
The appropriate difference equation for row relaxation in uniform mesh is given by 


ICC 
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and for column relaxation 


e*Ax L^‘i+lj ‘i-ljj 


(El 6) 


+^Vj!^ = o 

e hj 


Substituting equation (E4) into equation (El 6) and simplifying, yields the following 
equation for g 

fi- 2K(l - cos^i VAX'" + -^sin0, 

® , 2(l-cosg,) to^ 

^ Ay2 ■ ^ 


If M >w^ /ci then g <1 requires 


2co . « 2(1 -cos 02) 2K(l-cos0,) 

+ — T- sin0i <C 7 — o 1 1 — 5 

€ €*Ax * Ay- Ax^ 
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APPENDIX F 


SOLUTION PROCEDURE INVESTIGATION EXAMPLES 

F.l SOLUTION PARAMETERS 
F.l.l RELAXATION FACTORS 

The characteristics of the variation in solution convergence with relaxation factor magnitude 
are discussed in section 5.2.1. Some examples to illustrate these characteristics are included 
in this section. 

Generally, the coarser the mesh, the longer the initial convergence slope lasted, and thus the 
smaller the optimum overrelaxation factor (ORF), as shown in figure FI. Also, for a given 
frequency, as Mach number increased, the initial slope decreased, as shown in figure F2. Foi; 
a given Mach number, this slope decreased as the frequency increased. Generally, as the 
reduced frequency was increased, the ORF for most rapid convergence was reduced. As 
discussed later in this section, the size of the optimum ORF is also dependent on the 
relaxation procedure. 

Also shown in figure F2 are two examples of the effect on the convergence rate when the 
ORF is decreased during an analysis. The change is particularly noticeable for the M = 0.80 
curve when ORF is dropped to 1.4 from 1.85. 

An example of the effect of the magnitude of the underrelaxation factor (URF) is given in 
figure F3, and consists of a NACA 64A006 airfoil in harmonic pitch at a Mach number of 
0.9 and a reduced frequency of 0.06. The solution for the coarse, 25 x 20 grid converged 
relatively rapidly using an ORF = 1.6 and URF = 1.0. However, the solution for the finer 
42 X 30 grid diverged for URF = 1.0, but converged rapidly for URF=0.7. Variation of URF 
in similar calculations showed relatively little change in convergence rates as URF was varied 
between 0.7 and 0.3. 

As an alternative to the overrelaxation factors used so far, the pilot program was rewritten to 
accept complex ORF’s. Since there is little theoretical background for the use of complex 
relaxation factors, the procedure was to try several values to see what would happen. Here, 
only the real part of ORF was used in the calculation of ERROR. The results from using 
ORF’s of 1.6 ±0.1 i and 1.7 ±0.3i are shown in figures F4 and F5. Also included is the 
corresponding result for an all-real relaxation factor of 1.85, which appears to be optimum 
for this particular case. The ordinates for the figures are ERROR and the number of iterations. 
Although these calculations were not carried through to convergence, experience showed 
that the slope of this curve after it settled down (after some 30 or 40 iterations) was a valid 
indication of the rate convergence of the solution. Hence, an all-real ORF of 1 .85 appeared 
to be distinctly more efficient than any of the complex factors. ORF’s of 1.2 ±1 .2i were also 
tried, but ERROR diverged almost immediately. 

The alternative possibility of using different overrelaxation factors for the real and imaginary 
parts of the velocity potential was also examined. Figure F6 presents results for several 
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Figure FI. — Solution Convergence and Number of Points in Mesh 
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Flat plate with harmonically oscillating 
quarter-chord control 

CO = 0.06 

Mesh: 42 x 30; spacing factor: 1.4 






Figure F3. — Example of Effect of Underrelaxation Factor Variation on Solution Convergence 
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Number of iterations 


Figure F4. — Solution Convergence With Complex ORF's 
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Figure F6. — Solution Convergence With Separate ORF's for Rea! and imaginary Parts 
of the Velocity Potential 
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cases. T^e cases were run with ORF = 1.85 used for the real part and ORF = 1.3 and 1.0 for 
the unaginary part; then with ORF = 1.85 for the imaginary part and ORF 1.3 and 1.0 for 
the real, part. 

The mdst interesting result of this study is shown in flgure F7, where the convergence of the 
real and imaginary parts of the ERROR were plotted separately versus number of iterations. 
The periodicity as well as the phasing of the real and imaginary parts of ERROR was 
unexpected and seemed to be signiflcant. It perhaps helps to explain why convergence 
acceleration procedures such as those discussed in section 5.3.3. were unsuccessful. It would 
also seem to indicate that the complex ORF studies described above should be more successful. 

Finally, reference is made to appendixes C and E in which analyses were made of the relaxation 
solutions using time-like characteristics and the Von Neumann stability analysis. In both cases 
the addition of an imaginary term to ORF decreased the damping and thus did not aid in 
solution convergence. 

F.1 .2 GRID DISTRIBUTION AND SPACING 

In this section, examples are presented of pressure distributions in which there were variations 
in the spacing of adjacent finite difference intervals and in the number of points used to 
represent the flow field. The configuration used for these calculations was the flat plate with 
a Mach number of 0.8 and a reduced frequency of 0.06. As would be expected, in each case 
the representation of the singularities in the pressure distribution was improved by clustering 
the points about the singularities and by increasing the number of points. 

A program was developed for calculating the mesh-point locations for an airfoil with a control 
surface. The points were spaced according to the rule that adjacent intervals are X times the 
length of adjacent intervals with 2/3 < X <3/2. The user specified the scale factor X in each 
direction, the number of points over the wing surface (i.e., between x = -1.0 and x^), and the 
number of points in the vertical direction. In the vertical direction, points were placed 
symmetrically about y = 0 (the vertical dimension in a two-dimensional coordinate system) 
with no point at y = 0. Moving away from the y = 0 plane, each interval was X times the 
preceding interval, with the spacing set so as to put points on the upper and lower boundaries. 
The horizontal point spacing was more complicated, with points at the upstream and down- 
stream boundaries, and a point over the trailing edge, x = +1.0. Then points were equally 
spaced on either side of the leading edge, x = -1.0, and the control surface hingeline, x = x^. 

Also points over the wing (i.e., between x = -1.0 and x = x^) were symmetrically spaced about 
x_-l 

— jT If X’s of close to 3/2 are used, the finite difference points are clustered about the flow 
singularities at the leading edge and the hingeline in the flow direction and about the airfoil 
in the vertical direction. 

Examples of pressure distributions for several different mesh spacings are shown in figures F8 
through FI 1. Throughout the example, the area of mesh grid was held fixed and the total 
number of points nearly constant. The pressure distributions for the flat plate with an 
oscillating quarter-chord control surface are shown in figures F8 and F9. Corre^onding 
pressure distributions for the flat plate in pitch are shown in figures FIO and FI 1. In each 
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Flat plate with harmonically oscillating 
quarter-chord control surface 



Figure F7. — Convergence of Rea! and Imaginary Parts of ERROR 
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Figure F8. — Comparison of Solutions for Different Mesh-Point Spacings 
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Figure F9. — Comparison of Solutions for Different Mesh-Poin't Spacings 
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Figure F10. — Compar 


Flat plate in harmonic pitch 
M * 0.8, CO = 0.06 
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Figure F11. — Comparison of Soiutions for Different Mesh-Point Spacings 



case, the improvement in the representation of the singularity as the points were clustered 
about the singularities was very noticeable. 

A comparison of pressure distributions resulting when the number of mesh points was varied 
is shown in figures FI 2 and FI 3. The example is the flat plate with the oscillating quarter- 
chord control surface. The area of the mesh grid was constant and the points were placed 
according to the program briefly described earlier. These results are essentially directly 
comparable to the results presented in figures 4 and 5 of reference 1 for a much finer grid. 
There was a significant improvement in correlation with the linear solution as the number of 
points increased. It is felt that the 42 x 30 grid represented a good compromise between 
correlation with the linear result and overall economy, and thus would make a satisfactory 
basis for parameter variation calculations. 

F.1.3 EXTENT OF MESH 

An example of variations in the location of the upper and lower boundaries is presented next. 
This calculation was carried out at M = 0.8 and co = 0.06 for a flat plate with a harmonically 
oscillating quarter-chord control surface. The upper and lower boundaries (i.e., yj 
set at ±3.0, ±6.25, and ±10.0 in scaled coordinates. This is equivalent to ±8.99, ±18.5, and 
±29.6 in physical coordinates (i.e., 4, 9, and 15 chord lengths) for a 6% thick airfoil at 
M = 0.8. The number of mesh points was held constant for the three analyses using 42 points' 
flow-wise and 30 points in the crossflow direction. Figures F14 and FI 5 show the distribution 
of the jump in pressure coefficient across the section. The curves’for y, = ±6.25 and 
±10.0 are essentially the same except in the region of the leading-edge singularity. Here, a 
better representation of the singularity was obtained from the smaller yj value, apparently 
because the points in this case were more closely spaced about the wing section. The results 
at y 1 were slightly smaller in the amplitude for both the real and imaginary cases 

except in the region of the leading-edge singularity. Again, the representation of the leading- 
edge singularity benefited from the closer vertical spacing of the finite difference points. 

Figures FI 6 and FI 7 show plots of the velocity potential distribution for the variations in 
max discussed previously. The plots represent the variation in ip, (real and imaginary parts) 
in the vertical (crossflow) direction. The graph coordinates were set to emphasize the behavior 
of the velocity potential in the vicinity of the outer boundaries rather than the region next to 
the wing. The ipi for only two chordwise points are shown, but these curves are typical of the 
remaining curves. The curves for yj j^^x ±6.25 and ±10.0 lie very close together, while the 
corresponding curves for Vj of ±3.0 look quite different. Of course, it is the behavior of 
ipi adjacent to the wing that is most important, for this is where the pressure function is 
evaluated. However, for this case, the ipi distribution itself appears to have settled down by 
the time yi max “ ±6.25. 

F.1.4 SEQUENTIAL REFINEMENT 

As discussed in section 5.2.4, sequential refinement may be applied in terms of number of 
grid points, l^requency, or Mach number. Although all three forms were tried, it is sequential 
refinement with respect to mesh spacing that held promise of significant savings in terms of 
computer resources. Examples of sequential refinement with respect to both Mach number 
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Figure FI 2. — Comparison ofSoiutions with 
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Figure FI 3. — Comparison of Solutions with Difference in Number of Mesh Points 




Figure F14. — Comparison of Soiutions With Variation in Location of Upper and Lower Boundaries 
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Figure FI 5. — Comparison of Soiutions With Variation of Location of Upper and Lower Boundaries 
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Figure FI 6. — Variation of Velocity Potential With Location of Upper and Lower Boundaries 
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and frequency is shown in figure F 1 8. Depicted in this figure is the convergence for the case 
M = 0.9 and w = 0. 1 2 for three different runs. For the first, the initial velocity potential 
distribution was zeros; for the second, the converged solution for M = 0.85 and co = 0.12; 
and for the third, the converged solution for M = 0.9 and cu = 0.10. The difference in the 
number of iterations to convergence between the three runs was 25 out of 120, or about 22%, 
and is not significant unless the initial converged case is in existence. That is, it would not be 
worth generating the first solution in order to use it in the generation of the second case. 

The same is not true for sequential refinement with respect to the number of points while 
keeping the tptal area constant. Significant saving is not automatic, however, and several 
possibilities are sketched in the next two figures. The first figure (fig. FI 9) shows two possible 
convergence paths using a combination of a course grid with a fine grid. A straightforward 
linear interpolotation was used to obtain the initial values of for the iteration with the 
fine mesh from the last distribution of the course mesh. Note that there is a sharp peak in 
the convergence curves at the time the grids are switched. The difference between the two 
paths is the convergence value reached for the coarse mesh before the switch to the fine mesh 
is made. The computer cost is a function of the size and number of iterations performed with 
each mesh. For example, both paths reach the value C of ERROR in the same number of 
iterations. Path B is obviously more efficient than path A since a larger proportion of the 
iterations were made with the coarse grid. We found it advisable to carry the convergence of 
the coarse mesh below that of final convergence criterion. 

Another problem concerned the relative coarseness of the first grid to the second grid in a 
sequential refinement process. A second sketch (fig. F20) shows two convergence paths that 
differed as a result of different coarse grids. The evaluation of the total cost to convergence 
was more complicated than the preceding example. Here the efficiency must be evaluated by 
comparing the difference in cost of I 2 iterations of the fine mesh with the No. 2 coarse mesh, 
with difference in cost of Ij iterations of No. 2 coarse mesh with No. 1 coarse mesh. 

Thus it is very difficult to put quantitative numbers to the actual gain in computer efficiency 
using sequential refinement. We were satisfied with the results of using just two meshes to 
obtain final answers — a 25 x 20 mesh followed by a 42 x 30 mesh — for most of the examples 
of this report. 

Figure F21 presents the solution convergence for an airfoil using a 46 x 30 mesh. The upper 
curve shows the convergence starting with an initial (Pi distribution of zero. The lower curve 
shows the corresponding convergence starting with an initial distribution from a converged 
solution for a 1 7 x 20 mesh. The values at the points of the coarse mesh were linearly 
interpolated to provide values at the points of the fine mesh. The two meshes covered the 
same area. The point to be made is that the solution curves tend to merge for a large number 
of iterations and thus the advantage gained by using sequential refinement is lessened. As a 
side issue, the graph gives an indication of the effect of varying the scalar applied to the terms 
that are added to the finite difference equation to provide for convergence in the mixed-flow 
problem when using row relaxation. This factor, which is called CONPXT and is described in 
section F.2.1 , was set at 5 and 10 for this particular case. Also, the graph shows the effect of 
using an ORF of 1 .6 rather than 1 .85 for this problem. 
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Figure F18. — Sequential Refinement with Respect to Mach Number and Frequency 
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Figure F19. — Sequential Refinement With One Coarse Mesh and One Fine Mesh 
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Figure F20. — Sequential Refinement and the Effect of Coarseness of First Mesh 



Figure F21. — Examples of Effect of Starting Velocity Potential on Solution Convergence 
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F.2 SOLUTION PROCliSS 


F.2.1 ROW RELAXATION 


The finite dilTercnce equations presented by Eiders in reference 1 are written for column 
line relaxation. Rewriting these equations for row line relaxation consisted of simply putting 
the i+1, i-1, i-2 terms on the left-hand side of the equations and moving the j+I, j-l terms to 
the right-hand side of the equation. This change was accomplished in the program by minor 
DO loop changes and by rewriting the coefficient generator module. Although the mechanical 
aspects of tlie change were relatively straightforward, questions did arise regarding the effect 
on the solution cotivergcnce, a point made clear in discussions with Jameson. A general way 
of analyzing the iterative solution of finite difference ecpiations is discussed by Jameson in 
references 8 and 9. Itrielly the procedure is to treat the iterations steps as an artificial time 
coordinate. This permits rewriting the finite difference e(|uation as a differential equation 
including time, and the particular way the relaxation is performed is refiected in time- 
dependent terms. Thus a relaxation procedure may be considered as a finite difference 
equation for a time-dependent equation and the behavior of the corresponding iterative 
solution inferred from the resulting time-dependent equation. Both the subsonje (elliptic) 
and supersonic (hyperbolic) forms of the two-dimensional equation were analyzed; the 
details are included in appendix C. In summary, it was found that the elliptic form was 
convergent in its original form, whereas the hyperbolic form required the addition of two 
terms to achieve convergence. These terms are derived in the appendix and the following 
discussion will concentrate on the result of using row relaxation. 

Generally, row relaxation proved much more efficient than column relaxation. In addition, 
the row process was most efficient starting at the upper and lower boundaries and working 
in toward the wing surface, alternating the rows by taking one from the top section, one 
from bottom section, then one from the top section, etc. The one exception to this pattern 
was for combinations of Mach number and fretpiency for which convergence was marginal 
at best; then column relaxation proved superior to row relaxation. 

Figures F22 and F23 present examples of row and column solutions for comparison purposes. 
Figure F22 compares row and column solutions for a 42 x 30 grid for a fiat plate. It is 
believed that the solutions are carried out, initially at least, with a nearly optimum ORF. 
Figure F23 shows corresponding data for a 17 x 10 grid. In each case, row relaxation was 
significantly more efficient. 


Convergence in the airfoil case, i.e., for mixed flow, was obtained by adding the following 
terms to the finite difference equation. 
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Figure F22. — Convergence Comparison of Row and Coiumn Reiaxation Procedures 
for a 42 X 30 Mesh 





Figure F23. — Convergence Comparison of Row and Column Relaxation Procedures 
fora 17 X 10 Mesh 
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The values of the velocity potential with superscripts (n) are considered “new” values and 
are considered to be unknowns. The values with superscript (o ) are considered “old” values 
that were calculated in the preceding iteration and go on the right-hand side of equations 
(Al) and (A2). 

CONPXT and CONE6 are constants that may be varied from run to run. Column relaxation 
is convergent without these additional terms. A comparison between the solutions for row 
and column procedures (fig. F24) showed that row relaxation was again the most efficient 
process. The values used for CONPXT should be of the order 2 to 1 0 and for the value of 
COiNEb should be of the order of +0. 1 . As with ORF and URF, optimum values of CONPXT 
and CONE6 appeared to vary from case to case and, for maximum efficiency, should be 
determined separately for each calculation. 

There are several ways of running column relaxation for the two-dimensional case. The two 
most logical sequences appear to be: (1) starting from the upstream mesh boundary and, 
taking the columns in succession, moving to the downstream boundary; and (2) starting at 
the section trailing edge and moving forward, column by column, to the upstream boundary 
and then returning to the first column aft of the trailing edge and moving, column by column, 
to the downstream boundary. It is this latter process that worked best and is referred to as 
the “standard” column procedure. The former sequence will be referred to as simply “f>vd- 
aft.” An example using this latter solution sequence is shown in figure F25. From this 
particular example, it would appear that fwd-aft is more efficient than the standard sequence 
(see fig. F23). The fwd-aft sequence appears, in many cases, to be nearly as efficient as the 
standard column, but we found it generally to be less reliable. This was particularly true for 
combinations of Mach number and frequency that led to marginal solution stability. Also, 
the fwd-aft sequence utilized lower ORF’s, for several examples diverged when rerun with 
ORF’s that had been optimum for the standard sequence. It is interesting to note that working 
the solution sequence forward and aft from the the trailing edge does not seem logical for 
mixed-flow problems; however, in practice, it worked out very well. 

Similarly, there are other sequences for row solution besides the “out-in” sequence described 
above. These include, for example, “in-out” sequence (i.e., starting at the wing surface and 
working out, row by row, and alternating top and bottom), and just starting at one boundary 
(say the lower boundary) and solving successive rows to the upper boundary. A comparison 
of row in-out with row out-in is given in figure F26. 

F.2.2 DIRECT SOLUTION 

A version of the pilot program was also developed to provide direct solution for the interior 
velocity potential distribution. This concept was outlined by Ehlers (ref. 1). Generally, the 
iterative relaxation for the unknowns associated with the interior points was replaced by a 
solution of "the corhplete set of simultaneous equations at one time. There was still an 
iterative loop since the conditions on the outer boundaries of the mesh area were a function 
of the unknown velocity potentials adjacent to the airfoil surface. The modified pilot program 
was appropriate only for problems with subsonic steady-flow fields. 
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Figure F24. — Convergence Comparison of Row and Coiumn Relaxation for Mixed Flow 
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Figure F25. — Convergence Comparison of Row and Column Relaxation (Forward - Aft) 
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Figure F26. — Convergence Comparison of ‘ Un-Out" and "Out- In "Sequences for Row Relaxations 



The matrix equation for the direct solution may be written as 

i 

[AJ {X} = {R} (F2) 

f * . I 

where [A] is matrix of complex coefficients of order equal to the number of interior points, 
(X) is column matrix of unknown velocity potential, and {R} is a complex^collamn matrix 
introducing the conditions on the flow-field boundaries. If the number of niieslil points in the 
flow direction is i^^gx and the number in the crossflow direction is number 

of interior points and thus the order of equation (F2) is 

^ “ ^‘max ■ ^max ■ I 

The boundary values on the mesh region are dependent on the values of the vejpcity 
potential on the wing and wake; hence equation (F2) may be written as \ 

[A] {X} = {R(X)} '■ (F4) 

where [A] is a constant matrix and {R} is a complicated nonlinear vector valued 
function of (X} . The solution is obtained by iteration. That is, a {X^ is chosen and 
improved approximations to (X) are determined from 

[A] {X<"+J)} = {R(X")) , (F5)' 

I 

Note that for each fixed n, there is a linear equation of the form of equation (F2), and that 
as n changes, only the right-hand side changes, not [A] . 

The most efficient direct method for the solution of a linear system [A] (X} = {r 1 is well 
known to be the Gaussian elimination algorithm. The form of this algorithm that is, 
particularly suitable for several systems with the same coefficient matrix but with different 
right-hand sides is LU decomposition. In this scheme, [A] is decomposed, once and for all, 
into lower and upper triangular matrices [L] and [U], respectively, such that [L] [Ul = [A]. 
Then [A] (X) = {R} may be written as [L] [U] (X) = (R) and the solution found by 
solving, in turn, [L] (Y) = {R} and [U] {X} = {Y}. Note that since [L] is lower trialiigular, 
solution of [L] {Y} = {"R) iwolve^nly forward substitution; similarly, since [U] is upper 
triangular, solution of [U] {x} = {Y ) involves only backward substitution. Thus solution of 
{A] (x) = {r} is very fast once the decomposition [A] = [L] [U] is performed^ 

The main obstacle to solving large systems by direct methods is the large storage requirement. 
This is not due to the space required for [L] and [U] , since these may replace [A] as( they 
are calculated, but simply the storage required for [A] . The situation is improved soi^ie- 
what when [A] is a band matrix as in this case, but for a sufficiently fine mesh, the storage 
required may still exceed core capacity. To be more precise, when the coefficient ma,trix is 
banded, it may be shown that the band structure carries over to [L] and [U] . Thus it is not 
necessary to store all of [A] ; in particular, the O’s below the far subdiagonal and abpve the 
far superdiagonal need not be stored. (The O’s within the band must be stored since these 
locations are filled in the course of the decomposition.) Hence the total storage required for 
[A] is N [2(Jniax ‘ 1 1 • In addition to the storage for [A] we required N locatipns for 
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{R) and N locations for a scratch array used in the course of the decomposition. In sum, if 
the values were all real the total additional storage requirement, T, would be 


T = N* [2Gmax-2)+ n +2N (F6) 

= [20'max-2) + 3], (F7) 


but since the values are complex, we have 

T = 2N* [2Gmax-2) + 3]. (F8) 

We proceed to put (F7) in slightly different form, in which the dependence on the number of 
interior points is more transparent. 


Let 

Then 

which leads to 


^ ^'max ' 2)/(irnax ■ 2)' 


max 


^max ■ 2)^ ^^*max ' 2) Cimax ' 2) CN 


(F9) 


T = 4^^ • + 6N (FIO) 

Hence, we see that the storage requirement grows essentially as the number of interior points 
grows to the 3/2 power, for a fixed ratio of mesh spacings. For a mesh of 27 x 18 (i.e., with 
imax ~ 27 and j^ax “ 1 8), T = 28 000 words, which is within the core capacity of the 
CDC 6600. For a mesh of 42 x 32, T = 151 200 words, which is beyond the core capacity of 
the CDC 6600. 

F.2.3 CONVERGENCE ACCELERATION METHODS 

The relatively regular, uniform, and monotonic behavior of the pressure difference distributions 
with successive iterations suggested the use of convergence acceleration techniques. Also, these 
procedures were successfully applied in limited examples of steady transonic flows by Hafez 
and Cheng (ref. 13), and Martin and Lomax (ref. 14). Further, our studies showed relatively 
good behavior of both the unsteady velocity potential distribution and the ERROR with 
successive iterations. 

A rather straightforward program system was set up whereby three successive velocity 
potential distributions (or pressure distributions) could be calculated and saved. These would 
be for the n-2m, n-m, and n iterations, with both n and m being set by the user, A separate 
program is then used to generate a new distribution, which may be used as a starting point 
for a new sequence of iterations. 

The sample problem used for all the calculations was that of a two-dimensional flat plate with 
a harmonically oscillating control surface at M = 0.8 and a reduced frequency of 0.06, A 
mesh with 42 points in the streamwise direction and 30 points perpendicular to the flow was 
used. It is estimated that convergence for this case was reached in less than 400 iterations. 


142 



If 


The first convergence acceleration procedure used was an application "bf Aitken-Shanks 
nonlinear transformation (6^ process) to the velocity potential distribution (ref. 21). 


(n)*^ ^»ij 

‘ij .» 




>U 


(n-2m) 


iwrr 






(n-2m) 


(Fll) 


y 


The results were not encouraging. A typical result is shown in figure F27 and presents 
ERROR versus iteration. Equation (Fll) was applied using the velocity potential 
distributions from iterations 25, 50, and 75. The resulting velocity potential distribution 
was used as a starting point for 25 more iterations. The resulting ERROR’S are shown by 
the triangle symbols. After an initial sharp perturbation, the ERROR fell back to the level 
at which equation (Fll) was applied, and the convergence continued as if nothing had 
happened. 

As a second step, equation (Fll) was used on the pressure distribution rather than velocity 
potential distribution. Since the functional iteration solution process cannot be restarted 
with the predicted pressure distribution, successive results were monitored to see if they 
would converge. The process was to use increments of 25 (i.e., m = 25 in eq. Fll), starting 
with the 25th iteration. Thus pressure distributions were obtained using equation (Fll) with 
the results from iterations (25, 50, 75), (50, 75, 100), etc. Here again results were discouraging. 
The real part of the pressure behaved very well; the imaginary part had obvious problems up 
through the (100, 125, 150) set. Typical results are presented in figures F28 and F29. 


Examination of both the velocity potential and the pressures as calculated with equation 
(Fll) showed that the results were very sensitive to the input values when the values tended 
to lie along a straight line. Under these circumstances, for example, three positive values, 
obviously leading to more positive values, may well result in a negative value when used in 
equation (Fll). To eliminate this problem, the equation was rewritten in the following form: 




+ 


1 








- ^ 


(n-m) 





(F12) 


Then, a lower limit was set on the value of the denominator in the second term to restrain 
the resulting extrapolation within certain limits. An example of the use of equation (FI 2) is 
shown in figure F30. The results show no improvement over using just straight relaxation. 
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© Standard solution, ORF = 1.85 

V 75th velocity potential distribution 
calculated with equation (F11), 
n = 75, m = 25 


Flat plate with harmonically oscillating 
quarter-chord control surface 
M = 0.8, CO = 0.06 
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Figure F27. — Example of Solution Convergence Using the Aitken-Shanks 
Nonlinear Transformation 
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Figure F29. — Flat Plate With Harmonically Oscillating Quarter-Chord 
Control Surface — M = 0.8, oj = 0.06 
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Figure F30. — Fiat Plate With Harmonically Oscillating Quarter-Chord 
Control Surface -M- 0,8, cu = 0.06 
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